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ABSTRACT 


A comprehensive study for the inviscid numerical calculation 
of the hypersonic flow past a class of elliptic-cone derived 
waveriders is presented. The theoretical background associated 
with hypersonic small-disturbance theory (HSDT) is reviewed. 
Several approximation formulas for the waverider compression 
surface are established. A CFD algorithm due to Lawrence is used 
to calculate flow fields for the on-design case and a variety of 
off-design cases. The results are compared with HSDT, experiment, 
and other available CFD results . For the waverider shape used in 
previous investigations, the bow shock for the on-design condition 
stands off from the leading-edge tip of the waverider. It was 
found that this occurs because the tip was too thick according to 
the approximating shape formula that was used to describe the 
compression surface. When this was corrected, the bow shock became 
closer to attached as it should be. At Mach numbers greater than 
the design condition, a lambda-shock configuration develops near 
the tip of the compression surface. At negative angles of attack, 
other complicated shock patterns occur near the leading-edge tip. 
These heretofore unknown flow patterns show the power and utility 
of CFD for investigating novel hypersonic configurations such as 
waveriders . 
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Chapter I 

INTRODUCTION 


The design of trans-atmospheric and aero-space vehicles has been of great in- 
terest in recent years. The National Aero-Space Plane (NASP) is an example of the 
current effort to develop technologies to design a vehicle that will fly into orbit while 
taking off and landing like a conventional airplane. A simple generic diagram of an 
aero-space plane configuration is shown in Fig. (1-1). The generic configuration is 
divided into three parts: a forebody, a scramjet propulsion unit, and an afterbody. 
All of these parts are to be blended together as a smooth, interacting entity. Besides 
these basic parts, of course, there may also be wings and tails taking part in the 
overall configuration. 

One concept for the design of the forebody part of such a configuration is 
that of a waverider. A waverider shape offers a high lift and a small drag, and, in 
addition, provides favorable flow properties for the inlet of the scramjet propulsion 
unit. This investigation is directed towards a general study of the aerodynamics 
and flow fields associated with a special class of waverider configurations. 

The concept of a waverider has been around for some time. A good history of 
the original concepts and a discussion of some possible aerodynamic applications 
can be found in Kiichemann[l]. A waverider is constructed by identifying the stream 
surfaces of known supersonic flow fields as new solid surfaces that are connected in 
such a way as to form a new aerodynamic configuration. The flow field and aerody- 
namic properties of the waverider configuration are thus well known from the basic 
flows from which they were obtained. This basic flow and geometric configuration 
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are called the on-design conditions. When the waverider shape is held fixed, and 
either or both of the Mach number and orientation of the oncoming flow axe varied, 
the varied conditions axe said to be the off-design conditions. Whereas the under- 
lying concept of waveriders is that the on-design conditions are easy to calculate, 
the off-design conditions are usually very difficult to calculate. The development of 
Computational Fluid Dynamics (CFD) possibly provides the only practical means 
for studying the off-design properties of waverider configurations. 

The first CFD analyses of waverider flow fields were made by Jones et al. [2] 
and Jones[3,4]. These studies considered the elliptic-cone derived waveriders of 
Rasmussen[5] for which experimental results were available (Rasmussen et al. [6] 
and Jischke et al.[7]). The CFD calculations associated with Jones were simplified 
by using the full potential equations to describe the flow fields. Also, because the 
elliptic-cone waveriders are conical together with their inviscid flow fields even at 
off-design conditions, the similarity properties of conical flows could be fully uti- 
lized. In spite of the fact that the related elliptic-cone flow fields were slightly ro- 
tational, the irrotational potential-flow calculations provided fairly good agreement 
with experimental results and with the perturbation-analysis results upon which 
the elliptic-cone waveriders were based. This was probably due to the fact that 
the waveriders were actually not very slender, such that the viscous drag was small 
compared to the wave drag, and that the vorticity produced by small perturbations 
from axisymmetric flow was small. 

The purpose of the present study is to deal with the flow fields and aerodynam- 
ics of the elliptic-cone derived waveriders by CFD methods utilizing the complete 
Euler equations. Both the fluid dynamics and aerodynamics of the waverider flow 
fields are of interest together with the numerics of the CFD solution method. It was 
also desired to obtain results for viscous effects and heat transfer by utilizing the 
Parabolized Navier-Stokes (PNS) equations. Whereas the inviscid flow-field studies 
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were successful, the studies involving the viscous effects did not lead to success. In 
addition' to a broad description of the inviscid-flow results, some of the difficulties 
associated with the viscous-flow numerics will be addressed. 

Just recently, after the completion of the present work, several other papers 
have appeared that deal with the calculation of waverider flow fields by CFD meth- 
ods. Long[8] used the Euler equations to analyze an elliptic-cone waverider and 
a modified psuedo-waverider. Both configurations were fitted with afterbodies in 
their base regions. The only off-design results that were reported were for varied 
freestream Mach numbers. Jones and Dougherty[9] studied round-nosed, sharp- 
edged non-conical waveriders derived from axisymmetric conical flow fields. The 
Euler equations were used with emphasis on the grid generation to account prop- 
erly for the flow near the sharp leading edges. Only on-design conditions were 
treated. Liao et al.[10] studied elliptic-cone derived waveriders by a CFD Navier- 
Stokes simulation. These calculations considered the design Mach number fixed 
(Moo = 4), but varied the angle of attack. Although local skin-friction coefficients 
were calculated, the overall lift and drag were not. Because the shock is spread out 
over several grid elements, it was noticed that there was some flow spillage around 
the sharp edges of the waverider where in principle the shock should be attached. 
This was also noticed by Jones and Dougherty, as well as within the present investi- 
gation. Whereas these recent investigations have some regions of commonality with 
the present investigation, the present study is concerned primarily with the broad, 
overall description of the flow past a class of elliptic-cone derived waveriders, both 
on-design and off-design. The results are related to the numerics of the CFD cal- 
culation scheme, and thus this is of interest also. The intent and the results of the 
present investigation are different and more comprehensive than the aforementioned 
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The above three studies used existing codes based on unsteady 3-D equations 
with super computers like the Cray-2. On the other hand, in this study the steady- 
flow equations are utilized. Thus, the computational dimension can be reduced 
by one and the calculation can be made more effective and economical than the 
previous investigations. 

Some of the difficulties involved in this problem lie in handling the singularity 
produced by the sharp leading edge. The elliptic-cone waverider has two very sharp 
leading edges, and a bow shock is attached at the edges at the on-design condition 
in the ideal inviscid case. This can cause difficulties in both grid generating and nu- 
merical integrating procedures. For example, a slightly deteriorated grid structure 
near the tip will result in the divergence of the numerical calculations. Therefore, 
for this type of problem, both procedures should be carried out with great care in 
order to get successful results. 

Because of the particular boundary geometry and the importance of the flow 
near the body wall, the utilization of body-fitted coordinates is necessary. Among 
various methods for grid generation, an elliptic grid generator, which is a widely 
used scheme, is utilized, because it produces a very smooth grid. Although the 
overall geometry of the waverider looks simple, a lot of effort should be exerted 
to get a desirable grid structure especially near the tip. For that purpose Roberts’ 
stretchingfll], Sorenson’s[12] method, and Anderson’s adaptive grid method[13] are 
adopted. 

Numerous numerical algorithms and methods can be considered for numerical 
integration and discretizing the Euler equations which are used in this study. In 
this investigation, Lawrence’s STARS3D code based on his algorithm, which is a 
steady version of Roe’s[14] approximate Riemann solver and is in the class of up- 
wind schemes, is utilized to solve the hypersonic flows past elliptic-cone waveriders. 
The upwind schemes can be classified into two categories; Flux Vector Splitting 
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(FVS)[15] and Flux Difference Splitting (FDS). The FVS divides any flux into the 
positive and negative parts according to the sign of its relevant eigenvalues first and 
then discretizes them by using one-sided differences. The FDS determines a flux 
difference for two corresponding cell interfaces first and then discretizes the flux 
difference according to the relevant eigenvalues. Lawrence et al.[16] applied a Total 
Variation Diminishing (TVD) [17,18] schemes of Chakravarthy et al.[19,20] to the 
PNS equations to develop an algorithm which is in the category of the FDS. The 
TVD schemes have attracted much attention in the field of hypersonic flows, since 
they have some desirable properties associated with the handling of discontinuities 
like a shock in view of the stability and built-in dissipation. This algorithm uses 
the Finite Volume Method (FVM) to discretize the governing equations and an 
Alternating Directional Implicit method to integrate in space. The FVM, which 
is acquiring popularity recently over the conventional Finite Difference Method, is 
reported to have some advantages for the problems with irregular boundaries. In 
order to calculate inviscid flows, which were only possible for this problem with this 
code, the option for an inviscid flow was used. 

Despite the difficulties due to the singularities involved in the elliptic-cone 
waverider flow, the complete compressible inviscid solutions are sought without 
modifying the sharp geometry. To confirm the validity of numerical results for the 
hypersonic flow past waveriders, analytic approximate solutions obtained by means 
of the Hypersonic Small Disturbance Theory (HSDT) are calculated and compared 
with computational values. The HSDT result, which is possible to get only at the 
on-design condition, is used for checking the trends of the CFD results but not 
necessarily for seeking complete agreement. The main goal of this investigation is 
the study of the flow behavior around waveriders at off-design conditions where 
analytical solutions do not currently exist. 



Chapter II 

ELLIPTIC-CONE DERIVED WAVERIDER 

In this chapter we describe the development of an elliptic-cone waverider. Al- 
though such waveriders can be obtained from the basic supersonic flow past any 
elliptic cone, we shall be concerned with those cones having small eccentricities. 
In these cases, approximate results can be obtained by perturbation methods and 
hypersonic small-disturbance theory[5]. We shall outline the basic results here and 
leave the details for Appendix A. 

2.1 Perturbation Expansions 

We introduce spherical polar coordinates (r, 9 , <f > ) as shown in Fig. (2-1). A flow 
or a body is said to be conical if it does not depend on the radial coordinate r. 
We assume that a basic axisymmetric conical flow has been established. The basic 
circular-cone semi-vertex angle is 6 and the corresponding attached shock semi- 
vertex angle is (3. The uniform freestream flow, denoted by the subscript oo, is 
in the z direction, which is the axis of symmetry of the basic circular cone. We 
assume (along with Rasmussen[5]) that the basic axisymmetric flow is perturbed by 
introducing an elliptic-cone perturbation in the cone body shape. The elliptic-cone 
body is thus described by 

OM) = <5[1 - e 2 cos 2(f) + 0(e 2 )] 5 (2 - 1) 


6 
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where e 2 is a small-perturbation parameter that is a measure of the eccentricity. 
The corresponding shock shape is assumed to be 

6 s ((f>) = S[c — £292 cos 2(f> + 0(e|)]. (2 — 2) 

Here a = f}/8 is the ratio of the basic-shock angle to the basic-cone angle. The 
shock eccentricity factor g2 is to be determined as part of the perturbation solution. 

To account for a singular behavior in the solution near the surface of the per- 
turbed cone, Rasmussen [21, 22] introduced a new stretched variable 6 0 to replace 
the polar angle 0: 

Ojj-s s- <W) 

fi-6 9,(4,)- OM) ' 

Thus, for the new independent variable, 0q = 6 a.t the perturbed cone and 9 0 = (3 at 
the perturbed shock. The radial velocity u(8, cf>), polar velocity v(9, <f>), azimuthal 
velocity w(9, <f > ), pressure p(8 , tf>), density p(8, <f>), and entropy s(#, <f>) were then 


expanded in the series forms 

u(0, <f>) = u 0 (9 0 ) + e2U 2 (8 0 ) cos 2(f) + O(ej), (2 - 4a) 

v(9,(f > ) = v 0 (9 0 ) + e 2 v 2 (0o)cos2(f> + 0(e £), (2 - 4 b) 

iv(8,4>) = e 2 ia 2 (0o)sin2/ + O(el), (2 - 4c) 

p(6,(f>) =p 0 (^o)[l + e 2 P2(0 0 )cos2(f>\ + 0(4), (2 - Ad) 

p(Q,<t>) = po(^o)[l + e 2 R 2 (0 o )cos2/] + 0(e 2 ), (2 -4e) 

s(6,(f>) = sq + c v e2S2(9 0 )cos2(f> + 0(el). (2-4/) 


Here the subscript 0 denotes the unperturbed basic-cone quantity, and the subscript 
2 denotes the perturbation quantities. These expansions are a variation of the 
original expansions by Rasmussen and Lee[23]. An approximate solution for the 
above variables was obtained within the framework of hypersonic small-disturbance 
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theory. The results axe valid for large freestream Mach numbers and small cone 
angles 8 such that the similarity parameter Kg = M^S is held fixed at an arbitrary 
value. When Ks —* 0, the linear-theory limit is obtained, and when K& — > oo the 
hypersonic limit is obtained. 


2.2 Conical Stream Surfaces 

A streamline is described by the vector equation 

V x ds = 0, 


(2 - 5a) 


where ds is a differential element of distance measured along the streamline. In 
spherical coordinates this becomes 

dr rdd r sin 9d<f> 


u 


iv 


(2 - 5 b) 


Consistent with the first-order perturbations and hypersonic small-disturbance the- 
ory [21, 22], this reduces to 

dr = d9 0 = 6 0 dtf> (2 -5c) 

rV oo v 0 (6 0 ) e 2 w 2 (6 0 ) sin 2(f>' 

The second and third members of this equation describe the family of conical stream 
surfaces. The variables can be separated, and we have 

e 2 w 2 (# 0)^0 d(j> 


0ovo(0o) sin 2 <f>' 

An accurate approximation for vq(0q) is (see Appendix A): 


(2-6) 


Vo{0 0 ) = — Voo^O ^1 — • 


(2-7) 


If we denote by <f> = <f> a the azimuthal location where the conical stream surface 
intersects the shock do = /?, then Eq.(2-6) together with Eq.(2-7) can be integrated 
to give 


f 8 ° w 2 (9 0 ) d6 0 [* 

62 Jp 61-6* 4 


* d<t> 1 . 

. . , = -In 
sm 24 > 2 


tan <f> 
tan 4> a 


(2-8) 
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This is done in two different ways. The first way captures the behavior in the 
vortical layer, that is, near 9o = 9. It yields a simple result and is the basis for the 
actual construction of the waveriders in Refs. [6, 7, 26]. The azimuthal velocity 
tv 2 (9o) is a complicated function (see Appendix A), and the integral on the left 
side of Eq.(2-8) must be evaluated numerically as it stands. It is thus useful to find 
an approximation for this integral. Since the integral is singular at 9q = 6, most of 
the value of the integral comes when 9q is near 6, that is, near the surface of the 
perturbed cone. Thus we can set w 2 (9o ) = W2(S). This is a good approximation 
when W2(6 q) varies slowly across the shock layer; otherwise there may be some 
error in the shape near the leading edge. Correspondingly, in the integral we can 
set 9q — S 2 = (0 O + S)(d 0 — 6) = 2 6(9 0 — 6). When this is done the integral can be 
approximated simply, and we obtain 


( -e2W2(S)\ [#o - _ [ tan <f> ' 

V Wee J [in>\ “ 


(2-9) 


Let us define a new small parameter e by 


_ e 2 w 2 (S) 

e ~ • 


(2-10) 


This parameter is positive since w 2 is negative when e 2 is positive, and vice versa. 
Equation (2-9) can now be solved for 9 0 ((p) and written in the form 

a.fA\ r + — a 1 1 / e 


W) , , ^ [ tan<£ 

~r = 1+( ‘ T ~ 1) 


( 2 - 11 ) 


In terms of the physical variable 9(<f>), obtained from Eq.(2-3), this can be rewritten 
as 

9(<f>) __ 0 c (<f>) W) W) 1 tan^ l 1/e _ 

6 6 [ 6 6 J [tan^J ‘ 

Equation (2-12) describes the conical stream surfaces generated by hypersonic flows 
past elliptic cones. It is valid for small eccentricities, that is, for small t 2 . Since all 
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of the streamlines that comprise a conical stream surfaces pass through a common 
ray along the conical shock, they all have the same entropy. Thus, a conical stream 
surface is also a constant-entropy surface. 

It was first established by Ferri[24] that the surface of a conical solid body in 
a supersonic conical flow must be a constant-entropy surface. In general, except for 
axisymmetric shocks, the entropy must vary azimuthally immediately downstream 
of a conical shock. Thus, whereas Eqs.(2-4) may hold in the vicinity of the shock, 
they cannot hold in general near the surface of a conical body. Equations (2-4) are 
thus said to be the outer expansions, and another set of perturbation expansions 
must hold near the surface of the body. They are called the inner expansions, and 
they describe what is said to be the vortical layer. The vortical layer has been 
studied by means of matched asymptotic expansions [21,22,25]. It is found that the 
outer expansion for the pressure in Eq.(2-4d) is uniformly valid across the vortical 
layer to all orders, and the first-order azimuthal velocity in Eq.(2-4c) is uniformly 
valid to the lowest order. As might be expected, therefore, the approximation in 
Eq.(2-12) for a constant-entropy surface is also found to be uniformly valid across 
the whole shock layer. 

Although Eq.(2-12) is a uniformly- valid approximation for a conical stream 
surface across the shock layer, it is a composite solution, and thus it is not unique. 
There axe other approximate formulas that are correct to the same order of accuracy 
in e 2 . Another possible approximation is 


( 


e M y=i m ?+ 

0 0 


^ ^2 _ ( W) ^2 


tan <f> 


I 1 /' 


tcLXl j 


(2-13) 


This gives approximately the same result as Eq.(2-12) since (9+9 c ) and ( 9 a +9 C ) are 
nearly the same for slender cones in hypersonic flow. Equation (2-12) was used in 
Refs. [6,7] and Eq.(2-13) was used in Ref.[26]. A better representation for a conical 



11 


stream surface can be obtained by approximating the combination w 2 {6o)/{0q + <5) 
by a linear variation with 0 O from the body to the shock. Thus we would have 


W 2 (0q) W 2 (S) 


Voo(0o + s)~ 2^00 8 


+ 


w 2 (0) w 2 (8) 


Vooi/3 + 6) 2V 00 <5 J 


do 

J-6. 


(2 - 14) 


If we substitute this approximation into Eq.(2-8), then upon integration we obtain 


0o - 6 \ f 2 w 2 ((3) 

P — 6 ) CXP [\ ({T + 1)u12(<5) 



tan <j> 
tan ^ 3 


i/t 


(2-15) 


If the exponential term in this expression were ignored, then this result would be 
the same as Eqs.(2-9) or (2-11). Equation (2-15) is more accurate near the shock. 


2.3 Construction of an Elliptic-Cone Derived Waverider 

A conical constant-entropy stream surface, such as given by Eqs.(2-12), (2-13) 
or (2-15), is used as a lower compression surface in a new waverider configuration. 
The complementary upper freestream surfaces are taken to be a pair of triangular 
plane surfaces that pass through the axis of symmetry (z-axis) of the basic cone 
and intersect the elliptic-cone shock at the angle <j> = ±<f> 3 . A cross-section plane, 
perpendicular to the axis of symmetry, is shown in Fig.(2-2). The angle <f> 3 where 
the upper freestream planes intersect the lower conical stream surface is said to be 
the anhedral angle. A perspective lower rear view of an elliptic-cone waverider is 
shown in Fig.(2-3). 

For the on-design condition, the flow is aligned with the upper freestream 
surfaces, and the shock is attached to the sharp leading-edge tip at <f> = <f> 3 . This 
orientation is said to be at zero angle of attack a = 0. After the body shape is fixed 
based on the on-design conditions, the off-design conditions can occur by changing 
the freestream Mach number and/or by changing the orientation of the waverider 
(a 7 ^ 0). We shall not consider non-zero yawing (or sideslip) angle. Thus the 
flows considered will be right/left symmetric with respect to the vertical plane of 
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symmetry of the waverider. For Mach numbers less than the design Mach number, 
the bow shock will detach from the sharp leading edge of the waverider. 

2.4 Specific Waverider Shapes 

Two basic elliptic-cone waveriders are considered, denoted by A and B. For 
both cases, the on-design Mach number is M 0 0 = 4 and the eccentricity is e 2 = 0.1. 
The pertinent functions of Ks = M^S are determined from the formulas given in 

Appendix A. 

Table (2-1) 


Parameters 

Waverider A 

Waverider B 

6 

12° 

18.62° 

<f>3 

60° 

70° 

K s 

0.838 

1.30 

a 

1.62 

1.34 

92 

0.382 

0.597 

1/e 

6.132 

7.611 


The parameters for Waverider B correspond to the shape tested in Refs. [6,7]. 
Waverider A corresponds to the shape tested in Ref. [26], except that the eccentricity 
was e 2 = 0.05 in Ref.[26], which led to a corresponding value of 1/e = 12.264. This 
large value of 1/ e produces a very sharp leading edge on the waverider, and the initial 
efforts at calculating the flow led to numerical instabilities which were subsequently 
overcome. Consequently, the larger value of e 2 was selected for the present study. 

For each of the foregoing sets of parameters for models A and B, the compres- 
sion surfaces were calculated according to Eqs.(2-12) and (2-15). The models with 
compression surfaces calculated according to Eq.(2-12) are denoted by Al and B 1, 
whereas the models with compression surfaces denoted by Eq.(2-15) are denoted 
by A2 and B 2. The two compression surfaces for models Al and A2 are shown 
in Fig.(2-4). The two compression curves are considerably different in the concave 
region where the winglet blends into the body and outward towards the leading- 
edge tip. Figure (2-5) shows a comparison of the compression surfaces B 1 and B 2. 
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An analysis for the leading-edge tip angles is given in Appendix C. The waverider 
types A2 and B2 have considerably thinner leading-edge tips. The two types of rep- 
resentations will be used to test the influence of variations in the waverider shape 
on the flow-field properties. 



Chapter III 

GOVERNING EQUATIONS 


Three coordinate systems (Cartesian, Spherical, and Generalized coordinates) 
involved in the ensuing numerical investigation are shown in Fig. (3-1). The new 
Cartesian coordinates are labeled differently from those defined in Chaper II. The 
streamwise direction is x unlike z in the previous definition. For the numerical 
calculation we will use this new coordinate system throughout as well as the gen- 
eralized coordinate system. The rotated coordinates ( x',y',z ') are for aerodynamic 
force calculations which will be considered in Chapter IV. The waveriders studied 
here have conical geometry and thus to describe them spherical coordinates are also 
utilized. 

The Euler equations in Cartesian coordinates are introduced first in a strong 
conservation-law form. Although governing equations in this coordinate system 
are not utilized directly for numerical integration, their components axe involved 
in the computational process. For example, the initial conditions are imposed by 
means of Cartesian velocity components. Then the unsteady Euler equations are 
transformed in generalized curvilinear coordinates. This would be necessary for 
high-speed flow problems with a complicated boundary geometry. If the governing 
equations are expressed in a body-fitted coordinate system, then the imposition 
of boundary conditions becomes relatively easy, but the equations themselves have 
more complicated forms. In this investigation the steady Euler equations in the 
body-fitted coordinates are numerically integrated by means of a space marching 


14 



15 


shceme. For a steady problem this would be more economical and desirable with 
respect to the CPU time than using the unsteady equations. 

Finally, in the last two sections initial and boundary conditions for the wa- 
verider problem will be presented. Despite the unusual shape of the waverider it is 
easy to implemet those conditions. A computer code written based on a body-fitted 
coordinate system can be used for a large class of similar body shapes. That is, 
we can solve for elliptic-cone and waverider flows with the same code written for a 
circular cone based on the body-fitted coordinates. In other words, the imposition 
of the initial and boundary conditions for the flows of a circular cone, an elliptic 
cone, and a waverider, is the same inspite of the geometrical difference among them. 
Considering both of the body-fitted and Cartesian coordinates are body-fixed coor- 
dinates, we can distinguish the waverider flow problem from the other flow problem 
of such as a circular or elliptic cone by the generated grid, but not by the initial or 
boundary conditions. 


3.1 Governing Equations in Cartesian Coordinates 

The unsteady Euler equations without source terms (such as a body force) 
can be written in a strong conservation-law form in a Cartesian coordinate system 
(xi : x, y, z) as 


where 


f + Ef =0. (*:*.*«), 


»=1 


U = {p,pu,pv,pw,Et} , 

pui 


E' = { 


puiu + 6np 
pU{V + S i2 p 
PU{W + 8 iz p 

l (E t + p)ui ) 


(3-1) 
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E t = p e + ^(u 2 + v 2 + w 2 ) , 

where u,- are Cartesian velocity components ( u,v,w ), e is internal energy per unit 
mass, and denotes the Kronecker delta function. The ideal gas model is used 
here for simplicity: 

p=(-t- 1 )pe = — (3 - 2) 

All the above equations axe nondimensionalized by taking the freestream values as 
reference parameters. The nondimensional variables are, 



3.2 Governing Equations in Generalized Curvilinear Coordinates 

The Euler equations in a generalized curvilinear coordinate system can be ob- 
tained through the following transformation; 


C = £*(*» y, •?), (£* : r], C), i = 1,2, 3, 


where £(= £*) is taken as the radial direction. The t;(= £ 2 ) is in the crosswise 
direction and ((= £ 3 ) is in the normal direction to the waverider body wall. These 
coordinates will be generated numerically. The 77 and ( coordinates are orthogonal 
at the body boundary but not near the leading edges and for the rest of the flow 
field. Especially in case that the £ and ry coordinates are taken to be fixed at the 
body wall, the generalized coordinates' are said to be body-fitted coordinates. If we 
utilize the chain rule, the governing equations in a new coordinate system can be 
expressed in the strong conservation-law form again[ll]: 


dU_ dE { 

dt + ^ dV 
1=1 


= 0 , 




(3-5) 
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where 


where [.A] is defined by 



*4 


' E ' 



F 

ii 

«-h|- 

M 




. G , 


[g ) 


[A) = 


ir ty 6 

lx Vy V* 

Cr Cy C* 


and the Jacobian J is defined by the determinant of [.4] as 


(3-6) 


(3-7) 


d(x,y,z ) 


(3-8) 


If we examine Eq.(3-7), then we can see that all the information for the fluxes of 
Cartesian coordinates and the flows of generalized coordinates are exchanged by 
means of the metrics and the Jacobian. It is to be noted here that the metrics 
appear in a special form such that they are divided by the Jacobian. The metrics in 
this form, as well as the Jacobian, have geometrical meanings. The Jacobian is equal 
to the inverse of the cell volume and the metrics combined with the Jacobian are 
area element vectors, as will be explained later. Thus, once an appropriate grid is 
generated, both the metrics and the Jacobian can be determined from the geometric 
consideration of the grid. Therefore, the explicit form of the transformation function 
is not necessary. Even for a special case where we can identify a transformation 
function, it is not necessarily desirable to use the metrics and the Jacobian calculated 
by this function in a numerical analysis. 
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3.3 Steady Hyperbolic Equations 

In this study we use the steady Euler equations which axe obtained from Eq.(3- 
5). They axe written as 


3E dF dG_ 

d( + a; + ac, _ °' 


(3-9) 


To numerically integrate the governing equations in space they axe to be of a hy- 
perbolic type. The steady Euler equations become hyperbolic in the f direction, 
if we impose a restriction such that the eigenvalues of the pertinent flux Jacobian 
matrices ( dF/dE,dG/dE ) should be real[47]. The restriction requires supersonic 
inviscid flow in the streamwise direction everywhere. Fortunately, the waverider 
problem investigated here has such a characteristic except for the case with a very 
high angle of attack and thus we can use the steady equations with space march- 
ing technique. The ( coordinate in the steady equations plays the role of the time 
variable in the unsteady equations. 


3.4 Initial Conditions 

To start the numerical integration of the Euler equations by space marching, 
the initial flow values at a starting point of £ must be specified. Since the waverider 
has sharp leading edges and there is no inviscid subsonic region expected for a mild 
angle of attack, the supersonic freestream values can be used as initial conditions 
(ICs). Even though the final form of the equations is expressed in the body-fitted 
coordinates to fit the unusual geometry of the waverider, the imposition of ICs 
for the velocity can be achieved easily by the utilization of the Cartesian velocity 
components. The initial velocity components are: 

= Vqq cos a cos /?, 

vl i = Voo sin a, (3 — 10a, b, c) 

wl i = Voo cos a sin /?. 
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The effects of the angles of attack a and yaw /? axe implemented in this procedure 
by the above relation as shown in Fig. (3-2). For a special case without angles of 
attack and yaw where the freestream flow is in the x direction, the only nonvanishing 
velocity component is u\ i = Voo . Thus keeping and using the Cartesian velocity 
components in generalized curvilinear coordinates makes it easy to impose starting 
values, even for the particular shape of the waverider geometry. For the other flow 
variables such as pressure, density, and temperature, freestream values axe used. 
The initial converged solution is sought at the position x = 0.05 where the step- 
back procedure is taken until we meet any given criterion. Since we deal with the 
inviscid problem only and thus the flow is conical, we don’t need to integrate further 
downstream by space marching. 

3.5 Boundary Conditions 

The boundary conditions (BCs) without yaw (/? = 0) axe composed of three 
parts as shown in Fig. (3-3); the fax-field away from the bow shock, the upper and 
lower parts of the symmetry plane, and the body wall. For the case with angle of 
yaw, the symmetry condition cannot be used and thus the whole flow region must 
be included in the calculation instead of the half domain. 

Far-Field 

In a supersonic flow the flow downstream of a shock does not influence the 
upstream flow. Thus we can use the freestream values for the region which is far 
away from an expected bow shock located around the lower portion of a waverider. 
The outer boundaxy values (denoted by l = l ma x) are set equal to those of the 
freestream as 

^l mol = U frctstrcam- 


(3-11) 
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When the angle of attack is such that an expansion wave exists above the waverider, 
the far-field conditions must be imposed outside the Mach cone emanating from the 
apex of the waverider. 


Symmetry Planes 

The waverider configuration studied in this work is symmetric about the plane 
z = 0 and the flow around it will be also symmetric as long as its symmetry plane 
is aligned with the freestream flow direction ( fi — 0). To impose this condition two 
additional neighboring grid points are necessary across both the upper and lower 
symmetry planes, since the numerical algorithm used is the second order in the 
crosswise directions. Let 


(3 - 12) 


then 
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(3 - 13) 


Here the -1 in X plays the role of changing the sign of w across the symmetry line. 
The above relations show simply the reflection condition which is shown in Fig. (3-4). 

Wall 

Since the finite-volume method is used, the wall surface is composed of cell 
interfaces obtained by the primary grids. Thus, fluxes at the wall instead of flow 
variables should be specified. For inviscid flow the contravariant velocity W in the 
(■-direction is zero. That is, 


V ■ h = 0, =J> W = + %u; = 0, 

J 


where 


n = 


vc 

|vcf 


(3 - 14) 


(3 - 15) 
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As a consequence, all the inviscid normal fluxes to the wall vanish except for the 
pressure term. 



Chapter IV 

PROBLEM SOLVING PROCEDURE 


The first part of this chapter describes the problem concerning waveriders, and 
the solving procedure is presented. Then the characteristics and difficulties involved 
in the problem are set forth. Finally a simplified relation for aerodynamic forces is 
presented. 

4.1 Problem Description 

A waverider, a lifting body with high lift-drag ratio, is chosen as the forebody 
of a generic aero-space plane, as mentioned earlier. This study is to numerically 
solve the compressible flows past waveriders at supersonic speeds at both on- and 
off-design conditions. Because for supersonic flows the downstream flow does not 
affect the upstream flow, unlike subsonic flows, a waverider forebody can be studied 
independently of any downstream or afterbody configuration. 

Jones[2,3] solved numerically supersonic flows past waveriders by means of 
the full potential method[27] and obtained quite successful results. They are not 
strictly valid when the vorticity is sufficiently large, which occurs for a fixed finite 
eccentricity 62 when Ks becomes sufficiently large. In this investigation the complete 
Euler equations in a steady form are utilized. Since there are no restrictions on 
rotationality, unlike the potential method, we can calculate for arbitrary hypersonic 
flow regimes. 
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Composition of Program Packages 

The whole computer code is composed of three packages of program files. They 

are: 

(1) AGRID : generating grid 

(2) STARS 3D : solving for flows 

(3) CPLOT : plotting various graphs and obtaining approximate solutions for 
waveriders 

AGRID is the code that produces body-fitted coordinates. This can generate 
algebraic, hyperbolic, and elliptic grids. The details about elliptic grid generation 
are presented in Chapter V. The STARS3D code developed by Lawrence[28] at 
NASA, whose algorithm is described in Appendix B, is used to solve hypersonic 
flows past circular cones, elliptic cones, or waveriders with various conditions. This 
code was about 7,700 lines long, but it was reduced to 5,000 lines through modi- 
fication for grids and COMMON block structure and removing some unnecessary 
subroutines. The CPU time was about 2.25 hours on the IBM 3081 for an 83x41 
mesh. The CPLOT code is used for plotting the various variables obtained by 
running STARS3D through either a SURFACEII or a FORTRAN basic graphic 
routine. This also calculates approximate solutions by means of the HSDT. In 
addition PLOT3D developed by NASA is utilized for entropy contours and other 
plots. 

4.2 Problem Characteristics 

The numerical calculation of flows past waveriders involves difficulties that fall 
into at least three categories which are somewhat interrelated. The first category 
is associated with the singular nature of the waverider shape at the sharp leading 
edge. The second category is associated with the large gradients in the flow arising 
from the bow shock wave, the large gradients near the sharp leading edge especially 
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when the freestream Mach number is slightly less than the design Mach number, 
and possibly large gradients arising from vortical layers. The third category is 
the numerical algorithm itself, including the implementation of the wall boundary 
condition. 

The upper surface of the elliptic-cone derived waverider is a flat delta-shaped 
surface that intersects the lower curved compression surface along a straight line 
that forms the sharp leading edge of the waverider. The singular behavior of this 
leading edge presents difficulties in constructing a desired smoothly-varying grid 
structure near the leading edge. Several grid-construction schemes will be studied: 
an O-type grid, a fan-type grid, and an adaptive grid. If the grid is too skewed or 
if the grid cell volumes change size too rapidly, the numerical integration process 
will diverge. An adaptive grid could be considered to be used for the computation. 
This, however, would yield a grid skewness both for the sharp leading edge region 
and for the bow shock region. As a test case, an. adaptive grid was used with the 
hopes for an improved shock structure, but it did not lead to a successful converged 
solution. 

For the on-design condition, a conical bow shock is attached at the leading 
edge. The upper surface of the waverider is a freestream surface, and the lower 
surface is a high-pressure, high-density compression surface. As is the usual case 
for hypersonic flows, there is a large change in flow properties across the bow shock. 
For the waverider problem, however, there is an additional problem. This manifests 
itself dramatically when the freestream Mach number decreases slightly from the 
design Mach number. In this case the bow shock detaches from the leading edge, 
and a large gradient in the flow is established around the sharp leading edge of the 
waverider. The same sort of behavior would occur if the location of the leading 
edge were perturbed, either by design or by a slight error in the representation of 
the waverider surface. Since the numerical shock structure itself is spread over a 
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few grid elements, the representation of the flow near the sharp leading edge is very 
sensitive to small changes or errors. 

Since the asymmetric inviscid flow to be calculated is conical, a vortical layer 
exists that is pronounced near the symmetry plane of the compression surface. 
In principle, the compression surface of the waverider for the on-design flow is a 
constant-entropy surface since it is formed from a conical stream surface from a 
basic flow past an elliptic cone. There will thus be large gradients in the entropy, 
density, and temperature near the symmetry plane on the compression surface. 
There will also be corresponding large gradients in the entropy around the sharp 
leading edge. 

In this study we do not deal with viscous flows. If viscous flows were to be 
considered, however, the vortical-layer effects would be mitigated. In their place, 
nevertheless, large gradients due to viscous boundary layers would occur. These 
would produce further interactions of large gradients at the leading edge. 

The third category involves the numerical algorithm together with the way of 
the boundary condition imposition based on the finite-volume method. There is 
always some error involved in the numerical computation scheme. The boundary 
condition for the waverider surface pressure is imposed by means of the pressure 
at a half grid element away from the waverider wall. At the sharp leading edge, 
this is a potential source of difficulty. Except for the grid generation, the numerical 
algorithm used by the Lawrence code in this investigation is unmodified. For inviscid 
flows, pressure blips near the wall are detected even for the case of axisymmetric flow 
past a circular cone. For the related PNS viscous flows, wall pressure oscillations 
in the stream direction have been noted[28]. It is possible that a more accurate 
means of imposing the boundary conditions is needed. This will be discussed more 
later. The Lawrence code initializes the space marching operation by a step-back 
procedure using the freestream conditions as initial conditions. For the conical 
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inviscid problem this converges to the desired solution since the flow does not change 
in the radial direction. The situation is very much changed when viscous effects are 
involved since strong viscous interactions then occur at the vertex of the waverider. 

When the waverider is inclined to the freestream, that is, when it is at a 
nonzero angle of attack, it is possible that other flow features will arise that can 
cause numerical instabilities when the angle of attack is large enough. Such features 
might be the detachment of the bow shock, the occurrence of imbedded shocks, or 
the occurrence of vortical motions. 


4.3 Aerodynamic Forces 

The normal-force and axial-force coefficients (C/v, C A )axe defined as 

N „ A 


C N = 


C A = 


( 4 - 1 ) 


9oo*V Qoo Sb' 

where N and A are normal- and axial-forces respectively, S& is a waverider base 
plane area, and is the dynamic pressure of freestream flow. The base plane Sb 
is used for a reference area in order to compare with experimental data which use 
Sb also. 


Joo = jPooVi = j7PooA&, 


( 4 - 2 ) 


Sb = l 2 f tan 2 [6(<{>)]d<t>, 

J<t >= o 

where l is the length of the waverider. The pressure force acting on a waverider can 
be expressed by 

F = - / j s p dS = - / j(P~ Poo)fidS, ( 4 - 3 ) 

where n is an outward unit normal vector and S is a closed surface domain which 
embraces the whole waverider body. They are composed of three parts, i.e., a 
freestream upper-surface, a compression under-surface, and a base plane. Now 
consider the definitions of normal- and axial-forces, and metrics. 


N = j ■ F, A = i-F , dS = - VC, 


( 4 - 4 ) 
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where t and j are base vectors of body oriented coordinates. If the pressure of the 

base plane surface is assumed to have the freestream value, then the integration for 

__ 

the base plane vanishes by the second expression of F in Eq.(4-3). Therefore, the 
N and A in discretized form can now be expressed by 


^ = 2^[(Poo-p)y]i ) j, A = 2^[(p 00 -p)y] | -j, 1 < i < i'mar,! <j < jmax, 

(4-5) 

where the subscripts i,j are for tj,( coordinates respectively and j denote the 
wall pressures on the waverider upper and lower surfaces with excluding the base 
plane. Note that the signs of both and are negative for the lower compression 
surface where p > p 0 0 and thus N and A become positive. The Cn and C a for the 
whole waverider body now become 

^ _ 2 — P/Poo)^j]i,j 

Cn - 7 m^ s b 


(4-6) 


To get these coefficients the summation should be carried out for the whole forebody 
geometry. However, for an inviscid flow over a conical body like the waveriders 
studied here a simpler calculation can be done and its procedure is presented in the 
following. From the grid constructed for this problem, we can let 


AS, = J2 (V(/J)i (4 - 7) 

1=1 

For an inviscid flow we can apply the conical approximation such as 


Pi\,j ~ Pii,i — Pji 


(4-8) 



28 


where 1 < i i,t 2 < t mor . By means of Eqs.(4-7) and (4-8), the normal-force and 
axial-force coefficients can be obtained as 




Ca = 


‘OO 

2 


iM^Si 


£[i - 

Poo 


The ratio N/A can be easily obtained by Eq.(4-9) as 
AT/ A = ~ ~ 

S,'! 1 - (^)](^f 


(4-9) 


(4 - 10) 


This means that to calculate N/A for the forebody with conical geometry it is suf- 
ficient to consider a special portion of the forebody separated by two cross sections 
instead of the whole forebody. Cl and Cq for nonzero angle of attack are calculated 
by the following relation: 

cos a — sin a 
sin a cos a 



C "1 

< >, ( 4 - 11 ) 

{Ca J 


where Cn and Ca are based on the Cartesian coordinates (x,y, z) fixed to the 
waverider body. The rotated Cartesian coordinates for Cl and Cp, ( x',y',z ') are 
shown in Fig.(3-1). Now L/D can be easily obtained from Eq.(4-ll). 



Chapter V 

GEOMETRICAL CONSIDERATIONS 


A generalized curvilinear coordinate system is of common use nowadays in 
solving numerically fluid-flow problems with complicated boundaries. If we use 
body-fitted coordinates with two coordinates on the body, the boundary conditions 
can be imposed easily and accurately. Even though there exist grid generating codes 
such as GRAPE[12] and EAGLE[29], it would be more desirable to write a code 
for the waverider problem, because the elliptic-cone waverider studied here has a 
geometrical singularity and thus we have to control the grid as much as we need. 

Among various requirements of grid, smoothness and good grid control are ba- 
sically important. Thompson’s[30] elliptic grid generation, which can produce very 
smooth grids, is one of the most widely used schemes. The grids can be controlled by 
the inhomogeneous source terms called control functions in the Poisson’s equations 
which are used for grid construction. In determining control functions to generate 
a required grid for an elliptic-cone waverider two factors are taken into account: an 
orthogonal grid near the wall and a very sharp tip. To overcome the difficulty due 
to double singularities we introduce a Fan-type grid which was found quite efficient 
there. In order to capture a bow shock more accurately we may utilize an adaptive 
grid. To get a more desirable grid structure near the tip the boundary points axe 
redistributed by means of a stretching function. Finally, the metrics, Jacobian, and 
conical grid are described in brief. 
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5.1 Elliptic Grid Generation 

The set of equations for 2-D elliptic grid generation is 

V zz + Vyy — P(jh 0, 0* *b Cyy = Qixii 0* (5 — la, 6 ) 

subject to the Dirichlet BCs such as 

ri;z = /i(f7,Ci),y = yi(*?,Ci), 

(5 - 2 a, 6 ) 

r 2 ;z = / 2 (»?,C2),y = 02(77, C2), 

where /i,yi,/ 2 , and y 2 are functions specified for the inner (( = ) and outer 

(C = C2) boundaries. P and Q axe called control functions which control the grid 
shape. The solution of the above system has the form 

r? = r 7 (z,y), ( = ((z,y), (5-3 a, 6 ) 

in a physical domain. But our purpose is to use 77 , ( as new coordinates and accord- 
ingly they are to be used as independent variables instead of dependent variables. 
Therefore, we need the solution in the form 

2 = -z(77,C), y = y( 7 7> 0 , (5 -4a, b) 

in a computational domain. Figure (5-1) shows the physical and computational 
domains related to the elliptic grid generation. If we interchange the dependent 
and independent variables by means of the Jacobian theory (or chain rule), then 
the elliptic grid generator becomes 

az vv - 2 Pz v < + 7 z Cf = -^(Pz v -f Qz c ), (5 - 5a) 

ay^ - 2/3y nC + 7 y cc = --^(Py, + Qy <), (5 - 56) 

where 

Q = Z C +y<, P - + y^O 7 = *J + yJ- (5 — 6a, 6, c) 
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and the Jacobian J is defined as 

T _ dQ?.C) 

d( 2 , y y 


(5-7) 


5,2 Control Functions 

In determining control functions two factors are taken into account. One is 
grid control near the wall and the other is an adaptive grid. In this study the above 
two separated effects are combined together as 

P = Pw + Pa, Q = Qw + Qa, (5-8) 

where the subscript w stands for wall and a for adaptive. The Pw and Qw can 
be determined only by geometric constraints, while Pa and Qa are affected by flow 
solutions. 

Grid Control near Wall 

To control the grid near the wall two constraints by Sorenson[l2] are imposed. 
They axe orthogonality and the first grid spacing from the wall; 

d* y^yclw — 9, (5 — 9) 

i z l + y?) 1/2 it^ = Asi = specified, (5 - 10) 

where A 77 = A£ = 1 as usual and As* is the first grid spacing from the wall. Except 
for Ztf and y ^ at the wall all the other necessary derivatives can be determined if the 
inner wall boundary are specified. The second derivatives in £ are approximated by 
using the one-sided differences which contain previous iteration values. The source 
terms are assumed to have the following forms. 

PwiViO = p(v)^~ a< ', Qw{! hO = q(v)e- b<: . (5 -11a, b) 
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The above constraints can decide the p(r?), q(» 7 ) and the parameters, a and 6 , are 
specified. If a larger value of a or b is used, the effect of constraints decay very 
quickly. On the other hand for a smaller value it decays slowly. In solving the set 
of equations, the Successive Line Over- Relaxation (SLOR) method is utilized. Near 
the sharp corners, numerical instability may often occur. To remedy the problem, 
various methods including under relaxation and mixed finite difference schemes 
based on the sign of P and Q are used. 

Adaptive Grid 

To improve the resolution in regions where rapid flow variations occur, and/or 
to reduce the global error, an adaptive grid may be utilized. The basic idea of 
adaptive grid which provides automatic adjustment to the flow pattern can be 
obtained by equidistribution principle which for 1 -D case can be expressed by 


x n w = (A x)w = const. 


(5 - 12) 


This means simply that the mesh size is smaller when a weight function w is larger 
and vice versa. If we differentiate it once with respect to 77 and compare with 
the Poisson elliptic grid generating equation for 1 -D, it can be easily seen that the 
control function P is related to Thomas et al.[31] introduced control function 

as 

PA = (v 2 z + l 2 y)<t>(r],0, Qa = (C« +Cy)0(77,O- (5 -13a, b) 


Anderson[13] related these <f> and ip to the weight function 


w as 


1 dw _ 1 dw 
w drj ' xv dC, ’ 


(5-14 a, 6 ) 


In this study w is determined by 


<0 = 1 + A|Vp|/|Vp| max') 


(5-15) 
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where A is a constant and the number 1 is introduced to avoid infinite grid spacing 
in regions where the pressure gradient vanishes. The pressure gradient in Eq.(5-15) 
can be obtained by 


dp _ ( dp dp 
dz ^ y< dr] Vr} 
dp dp dp slr . 


(5 -16a, b) 


where y/g = j. Note that the pressure obtained by numerical integration is defined 
at cell center points and thus it is necessary to express it at primary grid points 
before calculating pressure gradients. If we use the weight function based on the 
above definition directly, then the grid might become rough. Thus, we adapt the 
following smoothing 

w(k, l ) = -^{4u;(fc, 7) + 2 w(k, l + l) + 2w(k + 1, 7) + 2 w(k — 1,7) + 2 w(k, 7 — 1) 

16 

+ w(k — 1,1 — 1) + w(k — 1, / + 1) + w(k + 1, / — 1) + w(k + 1,7 + 1)}. 

( 5 - 17 ) 

Since the source terms are defined by Eqs. (5-1 1,13), the Poisson equations are 
set up as 


ar„ - 2/3r, c + 7 r cc 


--{ 


JJ p(v)e ° c + 0 ^( 77 , C) 


ri, + 


— g( 7 )e 6C + 7^(77, C) 


r C 


}■ 


(5 -IS) 


where r = ( 2 , y) T . The Eq.(5-18) is the final form which is numerically solved to 
get a desired elliptic grid for a waverider. 


5.3 Boundary Conditions 

To generate body-fitted coordinates by solving elliptic partial differential equa- 
tions two boundary conditions at inner and outer boundaries should be specified. 


Inner Boundary : Waverider 

The compression under-surface of the elliptic-cone waverider, 9 W = 9 w (<f>), is 
constructed in spherical coordinates by either Eq.(2-12) or (2-13). 
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For — <f>s < <f> < (f> a , the lower compression surface is given by 

y = — x tan 8 W cos <f>, 2 = x tan 8 W sin <f>, 0 < <j) < <f> s , (5 — 19 a,b) 

y = —x tan 6 W cos <f>, 2 = —x tan 9 W sin <f>, —<f> 3 <<j><0. (5 — 20 a,b) 

For <f> 3 < <{> < 2 tt — <f> s , two planes which are aligned with the freestream flow are 
used to define a waverider configuration. 


y» < y < 0 , 2 


y, <<£<*•, 

y* 


Z s 

y 3 < y < 0 , 2 = — y, -7r < <{> < —<}>„. 
y» 


Outer Boundary : Ellipse 

The outer boundary is defined as an adjustable ellipse. 

(y + yo) 2 - 2 


6 2 


d — 2 ~ 
a 1 


where 


(5 -21a, b) 
(5 - 22a, b ) 


(5 - 23) 


y 0 = aiH, a = a. 2 H, b = azH, H = tan^|^ =0 . (5 — 24a, 6, c, d) 

where a l5 a 2 , and a 3 are parameters determining the outer boundary shape. 


Redistribution of Boundary Points 

The inner and outer boundary points around the tip are redistributed by using 
Roberts’ stretching[ll] which is more desirable than an exponential stretching, since 
for the region where the clustering is not wanted we can get more uniform stretching. 
The arc length, S , is defined by 




(5-25) 


where 0 < /? < 00 and ( = — 1— r. As /3 — * 0, more clustered grid near 5 = 0 can 

Sniaz 

be obtained. As (3 — ► 00 , they become uniform. Based on the arc lengths we can get 
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the clustered grid as we like. This is done by redistributing the x and y coordinates 
along the wall line by means of the interpolation according to the calculated arc 
lengths. 

Grid Types 

By changing the inner boundary values with slight modifications, we can obtain 
differrent types of grids such as O-type, Fan-type, and adaptive grids. 

The grid structure near the sharp tip plays a very important role especially in 
respect to the convergence in numerical calculations. The smoothness is the general 
requirement of any grid for body-fitted coordinates. To get such a grid the grid lines 
and the variation of grid cell volumes should be smooth. If we adapt O-type grid, 
which is shown in Figs. (5-2, 3) and obtained for the location x = 1 , the grid structure 
near the tip will be skewed so much that the grid line smoothness can be spoiled. 
This cannot be improved by clustering tj grid lines, while the cell volume variation 
is not too bad. Because of the convex surface with an infinite transverse curvature 
at the tip, constant £ lines near the tip region are clustered and thus grid spacings 
along 77 line which comes out of the tip are small. To adjust them we impose a 
positive source at the tip which is determined by trial and error. 

To improve the smoothness we introduce another type of grid, a Fan-type grid, 
shown in Figs. (5-4, 5) where several rays come out of the same tip point. The grid 
points are calculated at the position x = 1. This type of grid will not cause any 
problem in numerically integrating the governing partial differential equations, since 
we are using a finite-volume scheme and there is no flux into the body due to the 
zero cell area at the tip. However, it is desirable to get at least four triangular 
grid cells, considering that the numerical algorithm is 2 nd order accurate in the 
crosswise direction, and thus four cells are involved to calculate a flux. The grid 
used here has six triangular cells. This can represent the flow differences near the tip 
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properly. But there is a difficulty still in constructing this type grid, since increasing 
ray number will make the cell areas smaller and thus the smooth variation of ceil 
areas deteriorates. This can be alleviated by introducing a point source at the tip 
and adjusting the rj grid line spacings near the tip. The Fan-type grid is more 
appropriate to get a finer grid in the r] direction, while the O-type grid is more 
desirable for a finer grid in the £ direction. 

The adaptive grid for a Type- 1? waverider obtained by the Eq.(5-18) is shown 
in Fig.(5-6) and its magnified figure is shown in Fig.(5-7). Both figures are obtained 
for the location x = 0.05 where the numerical calculations are made. As can be 
seen, the grid is skewed very much near the clustered region. Even though the 
adaptive grid may improve the resolution of the shock in general, it doesn’t seem 
to be very desirable for the problem which has sharp leading edges where numerical 
instability can occur frequently. 


5.4 Jacobian and Metrics 

The Jacobian and metrics have their geometrical meanings and they can be 
determined by the grid geometry in Fig. (5-8). The indices (n, k, l) are for 
coordinates respectively. The vertices denoted by the solid circles (•) are called 
primary grid points and the cell center points denoted by X are called secondary 
grid points. 


Jacobian and Cell Volume Element 
Under the definition of the Jacobian J , 


7 = d(£,*7iC) 

" 5(x,y,z)’ 


(5-26) 


the J is related to the cell volume in the physical domain inversely, i.e., 


J = 


1 

AV’ 


(5 - 27) 
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where A£ = A tj = AC = 1 is set as usual. Therefore, the Jacobian can be obtained 
from the geometry of the grid. There can be various ways to determine the cell 
volume. In this investigation it is calculated by using the following relation[32], 


1 

J 


= i(AS”,r + A Sgl'l, + A s"« 0 2 ) . 0 X; /wl/2 - ^*-1/2, 1-1/2). 


(5 - 28) 


where r represents a position vector for the primary grid points shown in Fig. (5- 
8). For example, the point C in the figure is corresponding to /+1 / 2 * The 

magnitude of the area element A5£j' 1 is equivalent to the area of ABCD and its 
direction is in the normal to the area with the the positive f-direction. Whereas the 
fluxes are defined at the cell interfaces like A S^f 1 in the finite volume method, the 
flow variables are defined at the cell center point. In the Lawrence algorithm, the 
marching level is raised from (n + 1/2) to (n + 1) as shown in Fig. (5-8) . Therefore, 
the point (n + 1, k, l) is located at the center of the area element ABCD instead of 
the center point of the cell. 


Metrics and Cell Area Element 

The comparison of the discretized equation(B-ll) based on the finite volume 
method with its corresponding differential equation(3-5) in steady state can provide 
a useful relation which shows the physical explanation for metrics. For example, 
the £ metrics is expressed by means of the geometrical quantity, the area vector, as 

jWltf 1 = AS# 1 . (5 - 29) 


Therefore, once a grid are generated by any means, the metrics combined with the 
Jacobian can also be determined through their corresponding area element vectors. 
The area vector A5£'j’ 1 can be calculated by the cross product of two vectors as 


AS#' = \(AC x BD) 


- - r" +1 1 x tr n+1 

2 v *+l/2,J+l/2 7 Jt — 1/2,/— 1/2 ^ * Vfc — 1/2./+1/2 






_ f n + 1 'j 

k+ 1 /2,J — 1 /2 x ’ 


(5 - 30) 
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Note that the sign of the area element is decided by its corresponding outward unit 
normal n and the metrics in FVM are defined at cell interfaces. 

Conical Grid 

So far we described only 2-D grid, while the problem is 3-D. We can easily 
construct the desired 3-D grid through contraction or expansion of the generated 2- 
D grid, since the elliptic-cone waverider has a conical shape. Since the £ = constant 
plane is perpendicular to the x-axis, for the grid constructed here we have 


£„=£, = 0 . 


(5-31) 



Chapter VI 

CIRCULAR- AND ELLIPTIC-CONE FLOW SOLUTIONS 

In this chapter, in order to confirm the validity of the code STARS3D which 
is used for this study and also to check the accuracy of the HSDT approximate 
method, we calculate some inviscid flow variables for circular cones and compare 
them with known exact solutions. As a second means of comparison, the solution 
of the flow past an elliptic cone is obtained by numerically integrating the complete 
Euler equations. Considering any waverider investigated here is constructed from 
the known compressible flow field past an elliptic cone at a supersonic speed, it will 
be profitable to obtain a flow solution about the elliptic cone in order to facilitate 
the understanding of waverider physics. In studying hypersonic flows it is important 
to note that for any asymmetric supersonic or hypersonic conical flow there exists 
a vortical layer near the body surface where the entropy changes very rapidly. 
For checking the vortical layer we present entropy contours. Finally to enhance 
the understanding of the waverider flows especially near the leading edges, the 
shock locations by both the HSDT and numerical integration for elliptic cones are 
calculated and presented. 

6.1 Circular-Cone Flow 

As test cases we calculate the shock angles j3 and normalized wall pressure 
Pxxi/Poo of inviscid circular-cone flows for several values of the basic circular-cone 
half angle 6. For the calculation Moo = 4 and 7 = 1.4 are used. The numerical 
data are from STARS3D code and the approximate analytic data from the HSDT 
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by means of Eqs.(A-6,7,12). The numerical shock position is determined by the 
location of the largest pressure gradient. The shock is captured with one internal 
grid point and it is a very clean shock without any wiggles before and after the 
shock. For 8 = 12.5°, 17.5°, and 20.0° Sims tables[33] render exact solutions which 
are shown in Table (6-1). 


Table (6-1) 


8 

0 or p w 

Sims 

STARS3D 

HSDT 

AEs-e 

AEa-e 

12.5° 

0 

19.65° 

19.72° 

19.90° 

+0.36% 

+1.27% 


Pw! Poo 

2.307 

2.306 

2.331 

-0.04% 

+1.04% 

17.5° 

0 

24.08° 

24.16° 

24.07° 

+0.33% 

-0.04% 


Pw/Poo 

3.368 

3.365 

3.382 

-0.09% 

+0.42% 

20.0° 

0 

26.49° 

26.42° 

26.34° 

-0.25% 

-0.55% 


Pw / Poo 

4.006 

4.002 

4.014 

-0.10% 

+0.20% 


As can be seen, the maximum percentage errors of the numerical values to the 
exact values (denoted by AE^-e) for the surface pressure and the shock angles 
are 0.1% and 0.36% respectively for the given range of 8. These indicate that the 
computational results are very dependable, and thus that the STARS3D code can 
be utilized for other similar flow calculations with confidence. On the other hand, 
the maximum percentage errors of the approximate analytic values to the exact 
values (denoted by AEa-e) for both the surface pressure and the shock angle are 
slightly larger than 1%. 

Table (6-2) is for two values of 8 that are used to generate perturbed elliptic- 
cone flows. For the smaller 6 case the relative error of the analytic shock to the 
numerical shock (denoted by AEa-n ) is positive. But for the larger <5 case the error 
is negative. The wall pressure data of the HSDT show consistently larger values 
than the exact values in Table (6-1) and the numerical values in Table (6-2). 
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Table (6-2) 


8 

P or p w 

STARS3D 

HSDT 

AE A -n 

12.0° 

P 

19.30° 

19.52° 

+1.14% 


Pw/Poc 

2.217 

2.241 

+1.10% 

18.62° 

p 

25.19° 

25.07° 

-0.48% 


Pw/Poo 

3.641 

3.657 

+0.44% 


6.2 Flow Field past an Elliptic Cone 

Figures (6-1) show the solution to the compressible supersonic flow past an 
elliptic cone at Moo = 4 and a = 0°. It is obtained by numerically integrating the 
Euler equations. For the reason of comparison each physical variable is normalized 
by its corresponding freestream value. The elliptic-cone geometry is decided by the 
half basic cone angle 8 = 18.62° and small perturbation parameter = 0.1 (see 
Eq.(2-1)). This specific geometry with the given flow condition is introduced so that 
it is consistent with its corresponding waverider Type-B at the on-design condition. 
Figure (6-la) shows the grid (63x63 mesh) for the elliptic cone. Figure (6-lb) shows 
the disturbance velocity parallel to the y — z plane, that is, V = V v j + V z k. The 
dots indicate that the flow is in the x-direction without any disturbance due to the 
body. Between the dotted and arrow areas we can see a distinct boundary where 
a bow shock is located (roughly speaking). Through this shock the freestream flow 
is deflected. If we take a close look at the shock location, it can be seen that the 
shock stand-off is larger near the vertical minor axis than the horizontal major axis. 
This is in agreement with Eqs.(2-1,2) for the HSDT gi = 0.597 and a = 1.34. The 
azimuthal velocity component shown in Fig. (6-lc) goes to zero as <f> approaches 
0°, 90°, or 180°. Its maximum absolute value occurs at <f> = 45°. It is to be noted 
here that flow variables are not given at symmetry lines where primary grid points 
are defined. Since the FVM is used in this study, they are defined at secondary 







42 


grid points. The sign of w is negative as shown in Fig. (6-lc). That is because 
of the positive pressure gradient from the vertical minor axis toward the horizontal 
major axis, as can be seen in Fig. (6-ld). The pressure is a maximum at the major 
axis. Figure (6-le) shows the pressure distribution near the lower minor axis from 
the freestream side to the wall side. A very sharp and clean shock is captured with 
one internal grid point. Figure (6- If) shows a comparison of the pressure on the 
cone surface obtained by the numerical integration of the Euler equations and by the 
HSDT. At the minor- and major-axes the two sets of data are very close. The HSDT 
result in the middle of the figure is about 3 percent greater than the numerical result. 
Figure (6-lg) shows a comparison of the pressure distribution across the shock layer 
near the minor-axis symmetry plane. The HSDT shock is located more closely to 
the cone body than that by the numerical calculation. Figures (6-lh,i) show the 
pressure and Mach number contours. In both figures a bow shock is identified by 
coalesced lines. 

6.3 Entropy and Vortical Laver 

Figures (6-2a,b) show the entropy contours for the elliptic cones with 6 = 12° 
and 18.62°. Each of them corresponds to the Type- A or Type- 5 waverider. The 
constant entropy surfaces are getting closer as they approach the cone surface and 
they have a common tendency to embrace the body near the minor axis. This 
means the entropy near the wall changes very rapidly in the normal direction to 
the wall. A region with a rapid entropy change also has a large vorticity according 
to Crocco’s equation. The region of the rapid change near the wall is called the 
vortical layer. According to an analysis by Fern [24] vortical singularities exist at 
the minor axes of an elliptic cone where multiple values of entropy occur. In order 
to delineate the vortical layer distinctly, of course, we need accurate entropy values 
on the wall and symmetry lines. In this respect a finite-difference method (FDM) 
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would be more preferable for the purpose of capturing the vortical layer, since in a 
FDM both the wall and symmetry values are defined. A partially improved result 
for the vortical layer might be achieved by means of a finer grid in those regions. 

6.4 Shock Locations for Elliptic Cones 

Figure (6-3a) shows the shock locations for an elliptic cone (Type- A waverider 
generator) by the HSDT method and numerical calculation. It can be seen that the 
shock angle due to the approximate analytic method is greater than that obtained 
by the numerical integration of the Euler equations. Figure (6-3b) shows the shock 
locations obtained by the above. two ways for another elliptic cone (Type-5 wa- 
verider generator). For this case the HSDT shock is located inside of the numerical 
shock. These phenomena correspond to the results of Table (6-2). In the figures 
we can also see that the shock locations predicted by the HSDT and CFD methods 
are in reasonable agreement, the more so in the Type- A waverider generator. 




Chapter VII 

WAVERIDER FLOW SOLUTIONS 


In this chapter we present and discuss the numerical Euler solutions for super- 
sonic/hypersonic flows past elliptic-cone waveriders. The base on-design condition 
is for Moo — 4 and a = 0°. The chapter is divided into six parts. The first part 
discusses the on-design flows past the two A and B types of waveriders, and compar- 
isons are made with the results of the HSDT. The second part considers off-design 
flows where a = 0° is held fixed. The Mach numbers above the design condition 
are Moo = 4.5, 5.0, and 10.0. The Mach number below the on-design condition is 
Moo = 3. The third part considers off-design conditions for which Moo = 4 is held 
fixed and both positive and negative angles of attack are considered. The fourth 
part considers the off-design conditions for which both the angle of attack and Mach 
number are different from those of the on-design condition. The fifth part discusses 
the constant entropy surfaces for the on-design flows. The sixth part deals with the 
inviscid lift and drag of the waveriders as functions of freesteam Mach number and 
angle of attack. 

7.1 On-Design Flows : = 4, a = 0° 

We consider separately the on-design flows for the Type-A1,A2 waveriders in 
Figs. (7-1, 2) and then the flows for the Type-Bl, B2 waveriders in Figs. (7-3,4). Next, 
we compare the computational solution to the flow past the Type- .41 waverider with 
the analytic approximate solution based on the HSDT in Figs. (7-5) and to the flow 
past the Type-Bl waverider in Figs. (7-6). 


44 



45 


(1). Type-Al Waverider (K$ = 0.838) 

Figures (7-1) show various drawings for the numerical solution to the compress- 
ible flow past the Type-Al waverider at its nominal on-design condition. Figure (7- 

la) shows the O-type grid (83x41 mesh) that was used. The details of the typical 
O-type grids for the sharp tip region were presented in Figs. (5-2, 3). 

The disturbed velocity that is perpendicular to the x-axis is shown in Fig. (7- 

lb) . Its magnified flow field near the leading edge is shown in Fig.(7-lc). A bow 
shock is captured under the waverider, and it ranges from the symmetry plane to 
the leading edge region. Below the shock there are dots which indicate that the 
flow is perpendicular to the y — z plane and the freestream is not disturbed at 
all. In the region between the lower waverider wall and the dotted area we can 
see the velocity components in the y — z plane, which indicate the freestream flow 
is deflected through the shock. In the upper region above the waverider we can 
observe that the freestream is disturbed. The right side of this disturbed region 
is the extension of the bow shock from under the waverider, and its strength is 
weakened as it goes upwards. This weakened shock finally becomes a Mach cone in 
the upper symmetry plane. 

Figure (7-ld) shows a comparison of the bow shocks as calculated by means of 
the HSDT and as captured by the numerical integration. The computational shock 
position was defined by the locations where the largest pressure gradient along 
each constant tj line occurs. The shock captured by the computation is found to 
stand off from the leading edge instead of being attached, which should be expected 
for the ideal on-design condition. The attached shock can be expected only if the 
waverider geometry is based on the exact solution to the corresponding elliptic-cone 
flow and the solution to the flow past such a waverider is also exact. Thus, if either 
or both of the two conditions are not met, we cannot expect an attached shock 
in general. The waverider configuration used in this investigation is generated by 
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means of an HSDT approximation. For the flow past the elliptic cone that generates 
the waverider, we could already see the difference between the two shocks by the 
HSDT and the computation in Fig.(6-3a) of the previous chapter. Through the 
gap between the leading edge and the shock, the flow is appears to be spilling from 
the lower region to the upper region owing to the large pressure gradient in the 
circumferential direction. As a result, the flow in the upper region of the waverider 
is also disturbed. 

The normalized azimuthal velocity component w/Voo along the waverider wall 
is shown in Fig.(7-le). The horizontal axis denotes the ^-coordinate along the wall. 
The azimuthal angle <j> is measured from the lower symmetry plane anticlockwise. 
Obviously w goes to zero in the region at (f> = 0. As (f> increases, its magnitude with 
negative sign increases. In the middle of the figure we can observe that w increases 
from negative values to large positive values. In other words' the flow near the tip 
accelerates from the lower compression part to the upper freestream part of the flow 
field. 

Figure (7-lf) shows the normalized pressure distribution along the waverider 
surface p w /poo versus the horizontal axis z. For the computer computations, the 
z coordinate is measured at the location x = 0.05, where all the other numerical 
calculations are also carried out throughout this study. Since the flow and body 
are conical, the picture is similar in every vertical plane. The higher pressure line 
in the figure is for the waverider lower- compression surface, and the lower pressure 
line is for the upper-freestream surface. Strictly speaking, the pressure is for a half 
grid spacing above the wall. For the wall pressure the zero gradient assumption 
in the normal direction from the wall is used. As z increases, the pressure of the 
lower waverider surface also increases. This is mainly due to the larger deflection- 
angle effect from the minor axis to the major axis of the elliptic cone which is the 
waverider generator. For the same pressure line we can see the region where the 
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pressure decreases as z increases. This decrease is principally due to the pressure 
decreasing phenomenon of the elliptic-cone flow from the body to the shock as shown 
in Fig.(6-le). Near the leading edge, the pressure of the under-compression area 
increases as z increases, and it has the peak value right before reaching the edge. On 
the other hand, according to the analytical calculation which will be presented in the 
subsection Comparison with HSDT, the pressure for the same region decreases 
as z approaches the tip. To check whether the numerical pressure has the trend of 
the analytical calculation or not, we calculated the flow by using a more clustered 
grid near the tip and also Fan-type grids shown in Figs. (5-4, 5). But the pressure 
increasing trend remained almost the same. In other words, the behavior of the 
numerical solutions near the tip does not match that of the approximate solution, 
while the flow for the rest of the tip region remains nearly unchanged and shows good 
agreement with HSDT. This will be discussed in more detail subsequently. Figure 
(7-lg) shows the normalized pressure distribution near the lower symmetry plane 
z = 0. More precisely, it represents the pressure for the finite-volume elements 
whose left interfaces are located at the symmetry plane. The shock is found to 
be captured only with two internal grid points. It is a clean shock without any 
oscillations either after or before the shock. 

Figure (7-lh) shows the pressure contours. A bow shock is clearly seen in the 
lower compression part, and it stands off from the tip as stated earlier. The right 
upper area of the waverider is disturbed by both the bleeding due to the pressure 
jump near the tip and the extended shock. This disturbed region is expected to be 
extend upwards further as the angle of attack is increased. Therefore, for the case 
of high angle of attack, care must be taken so that the outer boundary is sufficiently 
far from the waverider to cover this disturbed region. If not, numerical instability 
cam occur. Figure (7-li) shows the Mach number contours which are similar to the 
pressure contours. 
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(2) . Type-A2 Waverider (K$ = 0.838) 

In order to determine the effects of a change in shape of a waverider on the flow 
fields, calculations were made for the waverider with shape Type- 42. This shape 
has a tip angle of only 1.44 deg, compared with 10.45 deg for the Type-41 waverider. 
The details for the tip angle axe described in Appendix C. Figures (7-2a) through 
(7-2g) pertain to the Type- .4.2 waverider, and these axe to be compared with the 
corresponding Figs. (7-la) through (7-lg) for the Type- 41 waverider. Figures (7- 
2a, a') show the Fan- type grid that is used for the computation. In spite of the 
sharpness of the leading edge, a very smooth grid is generated. The Cross-Plane 
velocity distribution is shown in Fig.(7-2b) and its magnified portion for the tip 
region is in Fig.(7-2c). We can see that the flow disturbance above the tip region 
is very small compared with the corresponding flow in Figs. (7-lb, c). The shocks 
by the HSDT and the numerical solution are depicted in Fig.(7-2d) and show an 
amazingly good agreement between them, unlike the case of the Type-41 waverider 
in Fig.(7-ld). From the comparison of the shock locations in Figs. (7-1, 2d) we can 
assert that the reason for the large discrepancy for the shock locations in the Fig.(7- 
ld) near the tip region is from the error of the waverider configuration at the sharp 
leading edge. The azimuthal velocity component is shown in Fig.(7-2e), and the 
pressure distributions along the wall and near the lower symmetry plane are shown 
in Figs.(7-2f,g). The waverider wall pressure distribution in Fig.(7-2f) shows no 
pressure peak at the tip, unlike the case in Fig. (7- If). The pressure distribution 
resembles the trend of the ideal on-design condition based on the HSDT. 

(3) . Type-Bl Waverider (K,s = 1.30) 

Figures (7-3) show various drawings for the numerical solution to the com- 
pressible flow past the Type-Bl waverider at its nominal on-design condition. The 
O-type grid (83x41 mesh) used for the waverider is shown in Fig.(7-3a). Since the 
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Type-51 waverider has a larger deflection 6 than the Type-Al, the flow of this type 
waverider has a larger hypersonic similarity parameter {Kg = M^S) and thus it has 
a more hypersonic feature. Figure (7-3b) shows the disturbed- velocity distribution, 
whose magnified portion for the tip region is shown in Fig.(7-3c). The disturbance 
region above the tip area can be found in the figure. This is located closer to the 
upper wall than that for the Type- ^41 shown in Fig.(7-lb). Thus, the possibility 
of interaction with the outer boundary is lessened. Due to this fact, for the case 
of high angle of attack this type waverider may give more stable results than the 
Type-Al. Since the Type-51 waverider is generated from the elliptic-cone flow 
whose numerical solution was presented in Chapter 6, ideally the flow pattern of 
the waverider lower compression portion should match that of the elliptic-cone flow 
pattern. We can see that the flow field in Fig.(7-3b) resembles the flow pattern in 

Fig.(6-lb). 

Figure (7-3d) shows two shocks from the approximate and numerical solutions. 
As the case in Fig.(7-ld), the numerical shock for this type waverider is located 
outside of the HSDT shock. Figures (7 — 3e ~ i ) show similar trends to those of 
the Type-Al waverider in Figs. (7 — le ~ t). 

For the flow past the Type-51 waverider, there are other data available for 
the wall pressure coefficient. Figure (7-3j) shows the comparison of the pressure 
coefficients C p by several researchers including the present study as functions of 
the azimuthal angle <f>. The experimental data by Jischke et al.[7] are plotted in 
the figure. There is a noticeable discrepancy between the experimental and present 
numerical data. The similar discrepancy between the experimental and the Euler 
numerical data by Liao et al.[10] can be observed. In fact, a complete agreement 
cannot be expected, since the experimental data are related to the real viscous 
flow. However, except for this basic difference of the viscous and inviscid flows, it 
is difficult at the present time to explain very well the reason for the discrepancy 
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in C p . In spite of this disagreement in magnitude between the experimental and 
computational data, the comparison shows very similar trends for the range of the 
given experimental data points. Furthermore, a very good agreement between the 
numerical results by Liao et al. and the present investigation can be found except 
for the very small region near the leading edge. The two results are obtained by 
completely different numerical integration methods. The above considerations on 
the trend and agreement lend confidence to the reliability of the present numerical 
results. In the figure the pressure coefficient by Jones[4] is also plotted. These data 
show better agreement with the experimental data for the lower values of but 
for the higher values they do not. The trend shows somewhat irregular variation at 
(f> > 30°. 

For the purpose of further checking the reliability of the numerical result we 
calculated the flow by using two different mesh numbers. The coarse mesh is 45 x 21 
and the fine mesh is 83 x 41. The comparisons for the different mesh numbers 
are shown in Figs.(7-3k,l). Figure (7-3k) shows the normalized pressure near the 
lower symmetry plane versus the normalized vertical length y/x tan S. The shock 
locations are almost identical. The wider shock structure denoted by the dotted line 
is due to the worse resolution by the coarse grid. Figure (7-31) shows the normalized 
waverider wall pressure distributions p w /poo versus the normalized horizontal length 
z/xtBJxS. As a whole, a good agreement except for the sharp leading edge region 
can be seen. The part of the reason for the discrepancy is from using the FVM, 
since the flow points for the grids of two different mesh sizes in the FVM cannot be 
identical. The other aspect is that the major numerical error occurs near the sharp 
leading edge. 
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(4). Type-B2 Waverider (K$ = 1.30) 

Figures (7-4a) through (7-4g) pertain to the type 5 2 waverider, and these are 
to be compared with the corresponding Figs.(7-3a) through (7-3g) for the type B 1 
waverider. The Type-52 waverider has a tip angle of 4.12 deg., compared with 
9.45 for the type B 1 waverider. Figure (7-4a) shows the Fern- type grid used for 
the computation. The cross-plane velocity distribution is shown in Fig.(7-4b), and 
the magnified region near the tip is shown in Fig.(7-4c). Although the shock is 
not attached at the tip, the stand-off distance is somewhat less than for waverider 
B 1. A comparison of the shock locations by the HSDT and the numerical solution 
is shown in Fig.(7-4d). The numerically calculated position is always somewhat 
outside the HSDT. This can be partly associated with the fact that, for this case, 
the conical shock angle at the tip is about 25°, and when the angles get this large the 
small-angle approximations of HSDT loose their accuracy. The azimuthal velocity 
component w is shown in Fig.(7-4e), and the pressure distributions near the wall 
and near the lower symmetry plane are shown in Figs.(7-4f,g). The magnitude of 
the azimuthal velocity is much less for this case, and the peak in the pressure at 
the tip of the compression surface has been removed. 

(5,6). Comparison with HSDT 

One way to confirm numerical results is to check them at a special condition 
where we know analytic solutions. The basic waverider-generating flow field past 
an elliptic cone can be obtained by a perturbation method. The details for the 
solution procedure axe described in Appendix A. Both approximate and numerical 
solutions for the Type-Al waverider are plotted in Figs.(7-5) and for the Type-51 
waverider in Figs. (7-6). The 0-type grids are used for each waverider. It should 
be kept in mind that the perturbation solution is an approximate solution, and 
the purpose of the comparison is to compare trends as well as magnitudes. As a 
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further reference, flow variables for unperturbed basic cones are also plotted; the 
polar angle 9 within the circular-cone shock layer is given by the same relation as 
the waverider wall equation 9 = 9{<f>) in Eq.(2-12). All the dependent variables 
are nondimensionalized by their corresponding freestream values and represented 
in spherical coordinates (r, 9, <f>). 

Figures (7-5) show the flow fields on the lower compression surface of the Type- 
Al waverider as functions of azimuthal angle The leading edge of the waverider 
is located at <f> = 60°. At <f> = 45°, the perturbed and unperturbed solutions are 
identical except for the azimuthal velocity component. 

The radial (r) velocity components of both circular and elliptic cones are shown 
in Fig.(7-5a). They have a common value at <f> = 45° as mentioned before. This 
is due to the perturbation term cos 2 <f> in Eq.(2-4a). At <f> = 0°, the r- velocity 
component of the elliptic cone is larger than that of the circular cone, as can be 
expected from the fact that the elliptic cone has a smaller deflection angle there 
than the corresponding circular cone. This trend remains the same for the region of 
0 < <j) < 45°. For the region 45° < <f> < 60° the opposite occurs by a similar reason. 
The trends of the analytic and numerical solutions are well matched. Even the 
magnitudes themselves are quite close to each other except for the tip region. For 
example, the numerical radial velocity component in Fig.(7-5a) near the symmetry 
line agrees with the analytic value within 0.15%. 

The polar (9) velocity components are shown in Fig.(7-5b). The difference 
between the components of the circular and elliptic cones is invisible on this scale. 
In other words, the 0-velocity component is insensitive to a slight disturbance on a 
circular cone. 

In Fig.(7-5c) the azimuthal ( <f > ) velocity components w are shown. The analytic 
value of w for the elliptic-cone flow has the term sin2<?i (see Eq.(2-4c)) which is 
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different from the rest of the dependent variables. The maximum value of the (f>- 
velocity component lies at 4> = 45° where the other perturbation variables vanish. 
Near the tip region, a relatively large discrepancy between the HSDT and the 
numerical calculation is detected. The much smaller numerical value of w from 
the computation is caused by the large pressure gradient in that region. 

The pressure distributions along the waverider wall are shown in Fig.(7-5d). 
At <j> = 0°, a larger pressure for the unperturbed cone is anticipated than for the 
corresponding elliptic cone. 

As presented in Appendix A, the pressure and azimuthal velocity component 
obtained from the outer perturbation expansions are valid for the entire shock layer 
including the inner wall regime. On the other hand, the rest of the flow variables are 
not strictly valid near the wall. This lack of validity of the outer solution is because 
of the existence of a vortical layer and a vortical singularity in the vicinity of the 
elliptic cone minor axis, which will be discussed in a later section. For a uniformly- 
valid solution we need to also use the inner expansions so that these phenomena 
might be taken into account. In line with this, it is interesting to note that for the 
symmetry plane region where the vortical layer and vortical singularity are most 
significant, the agreement of the approximate analytic and numerical values for the 
pressure and azimuthal velocity component in Figs.(7-5c,d) is very good. 

The density and temperature distributions according to the outer expansions 
are plotted in Figs.(7-5e,f). The agreement for those is not so good in comparison 
with the case of the pressure. 

Corresponding comparisons for the Type-51 waverider are made in Figs. (7- 
6). Figures (7-6a,b,c) show the normalized velocity components by the freestream 
velocity in the r, 9, and <{> directions. The comparisons for the radial and polar 
velocity components are very similar to those for the Type-Al waverider in Figs. (7- 
5a, b). The comparison for w in Fig.(7-6c) shows much better overall agreement 
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between the numerical and analytical values than that in Fig.(7-5c). Figure (7-6d) 
compares the waverider wall pressure. As for Figs.(7-5c,d), we can see the very 
good agreement for the pressure and the azimuthal velocity component in Figs.(7- 
6c, d) near the lower symmetry plane. In Fig.(7-6e) the outer-expansion density 
variations Eire plotted. A relatively large discrepancy between the numerical and 
the analytic values is observed. It would be helpful to recall that for the outer 
expansion solutions the pressure and the azimuthal velocity component Eire valid 
throughout the entire shock layer, but the density in Fig.(7-6e) and the temperature 
in Fig.(7-6f) axe not. 

7.2 Off-Design Flow : a = 0°. ^ 4 

In the following three sections, the various off-design solutions for the Type- 
BlO waverider are discussed. The O means that an O-Type grid is used. Even 
for the off-design conditions, as well as the on-design condition, the flow is coniced, 
since there is no characteristic length scale involved for inviscid supersonic flows. 
Therefore, for this problem it is enough to calculate the flow at a specific given 
position in the x-direction. In this study the waverider length is considered as 
unity and all the numerical calculations are carried out at x = 0.05. Thus, in the 
computation for conical flows savings in CPU time and storage axe significant. 

In this section the case for a lower Mach number than the on-design value is 
presented first and then the higher Mach number cases are presented. Finally their 
comparisons for the pressure distributions are made. 

(7). For Moo =3, a = 0° 

Figures (7-7) show various plots for the waverider flow at = 3 with no 
incidence. Figures (7-7a) shows the disturbance velocity parallel to the y — z plane. 
The shock stands off from the leading edge somewhat further than for the on-design 
case. In Fig.(7-7b) the normalized waverider wall pressure is shown as a function 



55 


of the horizontal Cartesian coordinate z. Near the tip a large pressure gradient is 
seen. Because of this large pressure gradient and the sharp comer of the waverider 
geometry at the leading edge, the flow at the upper surfce near the edge is disturbed 
significantly. For the flow at a Mach number below the on-design value, this kind 
of flow pattern near the tip area occurs quite typically. In Fig.(7-7c) the pressure 
across the shock layer at the lower symmetry plane is plotted. The normalized 
azimuthal velocity component w/Voo along the waverider wall is shown in Fig. (7- 
7d). The rapid change of w indicates the accelerated cross flow around the tip. 
In Fig.(7-7e,f) the pressure and Mach number contours are plotted and the shock 
stand off is clearly seen in those figures. 

(8) . For Moo = 4.5, a = 0° 

Figures (7-8) show various plots for the waverider flow field at Moo = 4.5 with 
no incidence. At this condition the bow shock is attached at the leading edge in 
Figs.(7-8a,b). The disturbance region is confined by the shock, which ends at the tip, 
to the under portion of the body. Figure (7-8c) shows the waverider wall pressure; no 
pressure peak occurs at the leading edge and no significant flow disturbance occurs 
on the upper freestream surface. This has a similar pressure trend to the idealized 
case obtained by HSDT. In Fig.(7-8d) the pressure near the lower symmetry plane 
is plotted, and the azimuthal velocity component along the waverider wall is shown 
in Fig.(7-8e). For the waverider shape B 1, this flow condition appears to be what 
should be called the on-design flow. 

(9) . For Moo = 5,a = 0° 

Figures (7-9) show various plots for the waverider flow field at = 5 with no 
incidence. In Fig.(7-9a) a bow shock is located inside of the tip. No flow disturbance 
is caught except for the lower portion of the waverider. The magnified flow for the 
leading edge is shown in Fig.(7-9b). In the figure we can see some other type of 
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flow disturbance from the bow shock to the tip. The main bow shock appears to 
terminate normal to the undersurface a short distance from the tip, and a very weak 
secondary shock emanates from the main shock, a short distance below the body 
surface, to the tip of the body. This shock interaction near the tip can be interpreted 
alternatively: For Mach numbers greater than for the on-design condition, the main 
conical bow shock will lie closer to the body. Except near the tip, the undersurface 
of the waverider is nearly the same as the originating elliptic cone, and the under 
portion of the bow shock corresponds to that for the elliptic cone. The flow near 
the tip, however, is governed by component of the freestream Mach number that is 
not parallel to the shock. At the tip 8 = 9 3 (<f> a ), this component is M± = M 0 0 sin 9 a 
and is parallel to the freestream surface. This component of the flow is termed 
the crossflow. The crossflow undergoes a deflection angle A at the tip, and this 
generates a weak oblique shock. The concave curvature near the tip causes a slight 
concave curvature in the oblique lip shock. The weak oblique lip shock and the 
strong conical bow shock intersect a short distance away from the lip. A strong 
Mach stem emanates from the line of intersection to a perpendicular location on 
the waverider body. The shock structure near the intersecting component has three 
legs which appear to have the shape of a “lambda”. We thus refer to the interacting 
shock as the lambda shock. We shall consider this further after examining the case 
for Moo = 10. A schematic diagram is shown in Fig.(7-10c). Figure (7-9c) shows 
the total pressure distribution for the tip region. Since the total pressure loss is 
also directly related to the entropy, this figure also indicates the entropy variation 
in the area. In Fig. (7-9d) the waverider wall pressure is plotted. There are two 
rapid pressure drops in the tip axea. The larger pressure gradient is caused by the 
main bow shock and the smaller pressure gradient is caused by the very weak shock 
located in the region from the main shock to the tip. The upper surface shows 
the undisturbed freestream pressure. In Fig.(7-9e) the pressure distribution across 
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the shock layer near the lower symmetry plane is drawn. The shock strength is 
stronger than that of the on-design condition as can be seen in Fig.(7-3g). Figure 
(7-9f) shows the azimuthal velocity component near the wall. Figures (7-9g,h) show 
the pressure and Mach number contours. The contours appearing to pass through 
the body surface near the tip in Fig.(7-9h) are due to the interpolation for the 
contour plots in the graphic package SURFACEII which is used for this figure. 
These also represent the lambda shock inside of the tip. 

(10). For Moo = 10, q = 0° 

Figures (7-10) show various plots for the waverider flow field at Moo = 10 with 
no incidence. Figure (7- 10a) shows the disturbed- velocity distribution on the plane 
perpendicular to the x-axis. The disturbance is limited to a small region below the 
waverider, and a bow shock is captured very close to the waverider compression 
surface. Near the tip region, a lambda shock pattern is established , and it is more 
pronounced than for the Moo = 5 case. It is shown more dramatically by the stag- 
nation pressure contours which axe shown in Fig.(7-10b) and which correspond to 
the Moo = 5 case shown in Fig.(7-9c). A schematic diagram of the lambda shock 
pattern is shown in Fig.(7-10c). The oblique shock from the tip is stronger for the 
Mqo = 10 case since the cross flow Mach number is larger and the cross flow under- 
goes the same deflection angle. Figure (7-10d) shows the wall pressure distribution. 
The rapid pressure jump on the lower surface occurs across the Mach-stem portion 
of the lambda shock that is normal to the surface. The upper surface has the undis- 
turbed freestream value. Figure (7-10e) shows the the pressure distribution across 
the shock layer near the lower symmetry plane. Even though the shock is quite 
strong, it is a very clean shock with only one internal grid point and no wiggles 
around it. In Fig.(7-10f) the azimuthal velocity component is depicted. 
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(11) . Comparisons for Flows at ^ 4, a = 0° 

Figures (7-11) show the comparisons of the pressure distributions at various 
frees tr earn Mach numbers = 3,4, 4.5, and 5 with the fixed angle of attack 
a = 0°. Figure (7-lla) compares the the pressure distribution across the shock 
layer near the lower symmetry plane. It can be seen that as M 0 e increases, the 
shock strength increases and the shock location moves toward the body. Figure 
(7-llb) compares the wall pressure distributions. The wall pressure for the upper 
waverider surface does not change so much except for the case of the Mach number 
below the on-design value. This implies that the upper freestream surface remains 
undisturbed for the higher Mach numbers above the on-design value. Especially at 
Moo = 4.5 the pressure variation near the tip resembles that of the HSDT prediction 
and the waverider upper surface has the freestream pressure. The wall pressure for 
the lower compression surface increases as the Mach number increases except for 
the leading edge region. 

7.3 Off-Design Flow : = 4. a ^ 0° 

In this section the cases for positive angle of attack are presented first and then 
the cases for negative angle of attack are presented. Finally their comparisons for 
the pressure distributions are made. 

(12) . For Moo = 4, a = +3° 

Figures (7-12) show various plots for the waverider flow field at = 4 and 
a = +3°. Figure (7-12a) shows the disturbance velocity in the y — z plane. Its 
magnified portion for the tip area is shown in Fig.(7-12b). In Fig.(7-12a) it can 
be seen that the flow region below the bow shock has uniform upward velocity 
which represents the freestream component due to the positive angle of attack. The 
bow shock stands off outside of the tip. The movement of the shock outwards 
from the on-design condition is caused by the increased effective deflection through 
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the positive angle of attack. Figure (7-12c) shows the pressure distribution along 
the waverider wall. The upper part of the figure represents the lower compression 
surface of the waverider, and the lower part represents the upper surface. We see 
that the pressure for the upper surface near the tip is smaller than the freestream 
value. The pressure drop is produced by an expansion wave around the sharp leading 
edge. We can detect a pressure rise near the symmetry plane for the upper waverider 
surface. This is produced by the formation of an imbedded weak conical shock that 
is needed to turn the flow that has expanded around the leading edge back to a 
direction parallel to the symmetry plane. This pressure rise was also observed by 
Jischke et al.[7] in their experiment for waveriders. The similar phenomenon can 
be seen by a supersonic flow past a delta wing at a positive angle of attack with 
a supersonic leading edge. Figure (7-12d) shows the pressure distribution near the 
lower symmetry plane. The shock is stronger than that of the on-design condition 
as in Fig.(7-3g). Figure (7-12e) shows the pressure distribution in the £ direction 
at the waverider top surface for the first grid cells whose left faces are at the upper 
symmetry plane. We can see the expansion that occurs from the freestream to the 
waverider upper wall. Fig.(7-12f) shows the azimuthal velocity component with the 
r] coordinate. 

(13). For Moo = 4, a = +10° 

Figures (7-13) show various plots for the waverider flow field at = 4 and 
a = +10°. As a whole, these figures present the characteristics of Figs. (7-12) 
more distinctly. Figure (7-13a) shows the disturbance velocity. The bow shock 
is located more outwards. The disturbance region extends more outwards than 
its counterpart in Fig. (7- 12a), and the undisturbed freestream exhibits now the 
uniform velocity components upwards clearly. The magnified disturbance velocity 
is plotted in Fig.(7-13b). The wall pressure is shown in Fig.(7-13c). We can see that 
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the cross flow around the leading edge affects the upper part of the waverider more 
forcibly. The pressure-rise region for the waverider upper surface is smaller than 
its counterpart in Fig.(7-12c). The shock near the lower symmetry plane and the 
expansion near the upper symmetry plane are shown in Figs.(7-13d,e). Both the 
expansion and shock become stronger than those in Figs.(7-12d,e). In Fig.(7-13f) 
the azimuthal velocity shows a positive peak value which indicates the accelerated 
cross flow at the leading edge. The azimuthal velocity w should be zero on the 
waverider upper surface. The w depicted in Fig.(7-13f) shows the values for the 
half grid spacing off from the wall. A small discontinuity is detected at rj = 75. 
This corresponds to the pressure-rise point near the upper waverider surface at 
z/xtan6 = 0.4 in Fig.(7-13c). 

(14) . For Moo = 4‘,a = —2° 

Figures (7-14) show various plots for the waverider flow field at Moo = 4 and 
a = —2°. At this condition the shock appears to be nearly attached at the tip 
as can be seen in Figs.(7-14a,b). In Fig.(7-14c) the pressure distribution along 
the lower wall does not show any peak value near the tip. In this respect the 
pressure trend for the lower compression surface resembles that of the HSDT result 
which represents the idealized on-design condition. It is interesting to note that 
at about this condition the maximum lift/drag ratio occurs, as will be presented 
in Section 7.6. Figures (7-14d,e) show the pressure distributions near the lower 
and upper symmetry planes respectively. In Fig.(7-14e) we can see the gradual 
compression from the freestream to the wall. Figure (7-14f) shows the azimuthal 
velocity component versus the rj coordinate. 

(15) . For Mqo = 4, a = —4° 

Figures (7-15) show various plots for the waverider flow field at M ^ = 4 and 
a = —4°. Figures (7-15a,b) show quite different flow patterns compared to their 
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counterparts in Figs.(7-12a,b) and Figs.(7-13a,b). First of all, the bow shock is 
located inside of the tip, which is caused by the reduced effective deflection angle 
through the negative angle of attack. Next, the upper disturbed region has the 
opposite flow direction to its counterpart in Figs. (7-12, 11a). That is, the upper 
disturbed velocity is in the direction of the leading edge. The gas near the leading 
edge flows downwards instead of upwards and the freestream has downward velocity 
components. In Fig.(7-15c) we can observe that, near the leading edge, the pressure 
on the upper surface is higher than the pressure on the lower surface. Thus, near 
the leading edge a negative lift is produced. The shock for the lower surface in 
Fig.(7-15d) becomes weaker than the corresponding situation with positive angles of 
attack in Figs. (7-12, lid). Unlike the cases with positive angles of attack, the upper 
surface along the symmetry plane experiences a compression from the freestream to 
the wall of the waverider instead of an expansion as can be seen in Fig.(7-15e). In 
Fig.(7-15f) the azimuthal velocity component shows negative values including the 
negative peak value near the leading edge. 

(16). For Mqo = 4, a = -8° 

Figures (7-16) show various plots for the waverider flow field at Moo = 4 and 
a = —8°. These represent the characteristics of the previous Figures (7-15) more 
dramatically. The bow shock in Fig.(7-16a) moves more inwards than that of Fig. (7- 
15a). Figure (7-16b) shows the downward flow more clearly. Figure (7-16c) shows 
the stagnation pressure contours where we can see the development of the upper and 
lower shocks. The shock below the waverider has a stronger strength near the tip 
area, and it becomes weaker as it approaches the lower symmetry plane. The shock 
above the waverider is nearly attached at the tip and weakens as it approaches the 
upper symmetry plane. In Fig.(7-16d) we can see the two areas delineated by the 
difference between the upper and lower wall pressures have approximately the same 
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size. This implies that the lift is alm ost zero. The situation of zero lift will be seen 
again in Section 7.6. The compression by the shock in Fig.(7-16e) becomes weaker 
than its counterpart in Fig.(7-15d), and the compression in Fig.(7-16f) becomes 
stronger than its counterpart in Figs.(7-15e). In Fig.(7-16g) the azimuthal velocity 
is shown. The azimuthal velocity component is always negative on the lower surface. 

(17) . For Moo = 4, a = -12° 

Figures (7-17) show various plots for the waverider flow field at Moo = 4 and 
a = —12°. Figures (7-17a,b) show somewhat irregular flow patterns for the wa- 
verider lower portion near the tip area. In Fig.(7-17c) we can observe that the 
lower surface pressure is smaller than that of the upper surface. As a result, a 
negative lift is generated, which will be presented again in Section 7.6. Figure 
(7-17d) shows the pressure distribution in the shock layer under the body near the 
symmetry plane; the flow downstream of the bow shock undergoes a gradual com- 
pression. Figure (7-17e) shows the pressure in the shock layer above the body near 
the symmetry plane. Figure (7-17f) shows the azimuthal velocity component which 
has a somewhat irregular variation. 

(18) . Comparisons for Flows at = 4, a ^ 0° 

Figures (7-18) show the comparisons of the pressure distributions at various 
angles of attack a = — 8°, — 4°, 0°, +3°, +10° with freestream Mach number fixed 
at Mqo = 4. Figure (7-18a) compares the pressure distribution and shock locations 
near the lower symmetry plane. It can be seen that as a increases, the shock 
strength increases and the shock moves toward the body. In this respect the effect 
of the increasing a is similar to that of the increasing Moo as can be seen in Fig. (7- 
11a). Figure (7-18b) compares the wall pressure distributions. For the negative 
angles of attack the pressure lines for the lower waverider surface intersect those 
for the upper waverider surface. The areas of the left hand side of the intersection 
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denote positive lift, but the areas of the right hand side denote negative lift. This 
is the reason for the decreasing lift at negative angles of attack. At the positive 
angles of attack the pressure near the upper symmetry plane becomes higher than 
that near the tip region, because a weak imbedded shock develops on the upper 
surface. We can also see the pressure near the lower symmetry plane increases in 
comparison with the tip region, as a increases. This is because the flow deflection 
near the symmetry plane due to the high angle of attack is greater than the flow 
deflection due to the perturbed elliptic cone near the major axis, as a increases. 

7.4 Off-Design Flow : ^ 0. a ^ 0° 

In this section the cases of different Mach number and angle of attack from 
the on-design values are presented and then their comparisons for the pressure 
distributions are made. 

(19). For Moo = 5, a = +4° 

Figures (7-19) show various plots for the waverider flow field at il/oo = 5 and 
a = +4°. Figures (7-19a,b) indicate that the bow shock is slightly detached from the 
tip. The larger freestream Mach number than the on-design value causes the shock 
to move inwards. On the other hand, the positive angle of attack causes the shock 
to move outwards. These two opposite effects work together to result in moving 
the shock nearer to the tip. In Fig.(7-19a) the freestream flow shows the upward 
velocity component. Figure (7- 19c) shows the wall pressure distribution. As shown 
in the cases of other positive angle of attack with the fixed Mach number Moo = 4, 
we can see a small pressure increase near the upper symmetry plane. Figure (7-19d) 
shows the pressure distribution and the shock near the lower symmetry plane. Due 
to the larger Moo than the on-design value and the positive angle of attack the shock 
becomes very strong. The pressure distribution near the upper symmetry plane is 
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shown in Fig.(7-19e), which shows an expansion from the shock to the body. The 
azimuthal velocity component near the body surface is shown in Fig.(7-19f). 

(20) . For Moo = 5, a = —4° 

Figures (7-20) show various plots for the waverider flow field at = 5 and 
a = —4°. In Figs.(7-20a,b) we can see that the shock is located inside of the tip. 
Both the larger Mach number and negative angle of attack cause the shock to move 
inward. The freestream portion shows the downward velocity component. For the 
region near the top ridge of the waverider we can observe that there is no velocity 
component in the y — z plane as in the case in Fig.(7-19a). Whereas the vertical 
flow is directed toward the ridge, the flow parallel to the waverider upper surface 
is directed outward from the ridge. The cross flow near the leading edge shows a 
smooth variation. The wall pressure line shown in Fig.(7-20c) intersects itself at the 
lower right side of the figure. From this we can see that the upper surface pressure 
near the leading edge is larger than that of its corresponding lower surface with 
the small pressure gradient. Therefore, the flow near the tip is downward with a 
smooth variation as seen in Fig.(7-20a,b). If we compare the pressure distribution 
in the lower symmetry plane in Fig.(7-20d) with that in Fig.(7-2g), it can be seen 
that there is no big change in the shock strength. That is because the larger Mach 
number causes an increase of the shock strength, whereas the negative angle of 
attack causes a decrease. The pressure distribution near the upper symmetry plane 
is shown in Fig.(7-20e), and the azimuthal velocity around the body surface is shown 
in Fig.(7-20f). 

(21) . For Mqo = 3, a = +4° 

Figures (7-21) show various plots for the waverider flow field at Moo = 3 and 
a = +4°. As can be seen in Fig.(7-21a) the disturbance area of the upper flow 
region extends far away. Both the larger angle of attack and smaller Mach number 
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than the on-design values cause the disturbance area to move more outwards. At 
this condition numerical instability can often occur due to the extended disturbance 
and the irregular flow pattern near the upper tip region as shown in Fig.(7-21b). 
In Fig.(7-21c) we can see a large pressure gradient near the tip. In the figure we 
see that the pressure change of the upper surface is significant, but the pressure 
variation for the lower surface is mild. In Fig.(7-21d) the pressure distribution for 
the lower symmetry plane is illustrated. For the pressure distribution near the 
upper symmetry plane, shown in Fig.(7-21e), the flow expansion field is fairly large. 
From the right upper portion of the figure, we can see that the pressure near the 
outer boundary seems to be influenced by the boundary which is not taken large 
enough to embrace the disturbance region. The azimuthal velocity component near 
the body surface is shown in Fig.(7-21f). A small discontinuity can be detected at 
7 / = 65, as happens usually for the cases with a > 0° so far. 

(22). For Moo = 3,a = -4° 

Figures (7-22) show various plots for the waverider flow field at = 3 and 
q = —4°. Figure (7-22a) shows the flow field perpendicular to the x-axis. For the 
upper part the weakened extended shock and the compression due to the negative 
angle of attack are combined together. In the figure we can see that the outer 
boundary near both the lower and upper symmetry planes does not embrace the 
whole disturbance region. In order to get more desirable solution we need to extend 
the outer boundary more outwards. In Fig.(7-22b) the magnified flow field for the 
sharp leading edge area is shown. The flow above the upper waverider surface, which 
is somewhat parallel to the wall, is directed toward the leading edge. Before reaching 
the tip the flow is deflected upwards a little bit due to the pressure gradient in the 
region. Figure (7-22c) shows the wall pressure distribution. We can see the pressure 
rising near the leading edge for the upper wall. In Figs.(7-22d,e) we see that both 
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the lower and upper surfaces near the symmetry planes undergo compressions. The 
compression for the upper surface is due to the negative angle of attack and possibly 
due to the extended bow shock. The other compression for the lower surface is due 
to the flow deflection caused by the waverider geometry, as is usual for the cases so 
fax. Figure (7-22f) shows the azimuthal velocity component near the body wall. 

(23). Comparisons for Flows at ^ 4, a 0° 

Figures (7-23) show the comparisons of the pressure distributions at the two 
angles of attack a = ±4° and the two freestream Mach numbers Moo = 3, 5. Figure 
(7-23a) compares the pressure distributions and the shock locations near the lower 
symmetry plane. As far as the shock strength and location are concerned, increasing 
the angle of attack and increasing the Mach number exert similar effects. Thus, at 
the conditon of Moo = 5 and a = +4° the shock strength is strongest and the 
shock is located nearest to the waverider body wall. For the case of Moo = 3 and 
a = —4° the opposite phenomenon happens. Figure (7-23b) compares the wall 
pressure distributions. The closed area in the figure by a pressure line is a measure 
for the lift. We can see that the lift at Moo = 5, a = +4 is very large but at 
Moo = 3, a = — 4 it is very small. 

7.5 Entropy Distribution 

Figures (7-24) show the constant entropy contours for the flows past the Type- 
5 1 waverider with the two different types of grids. Figure (7-24a) is for the Type- 
51 waverider with a Fan-type grid (denoted by BlF) and Fig.(7-24b) for the same 
waverider with an O-type grid (denoted by 510). The comparison of the three 
plots of Figs.(6-2b) and (7-24a,b) provides us one feature of entropy increase. The 
entropy for the waverider flow in Fig.(7-24a) with a Fan-type grid is increased in 
comparison with its corresponding elliptic-cone flow in Fig.(6-2a). This is mainly 
due to the sharp leading edge which makes the flow change very rapidly. Further 
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entropy increase can be seen in Fig.(7-24b) which is for the same waverider with an 
0-Type grid. The major difference between Figures (7-24a,b) is the grid structure 
near the tip. The 0-Type grid is more skewed than the Fan- Type grid near the tip. 
The entropy production caused by the skewed grid may be explained in terms of a 
numerical viscosity. For viscous flow the shock location moves outwards compared 
with that of the corresponding inviscid flow due to the increased deflection by the 
momentum displacement. Accordingly this numerical viscosity may cause the shock 
to move outwards more or less. Recall that for both the Type- A and B the shocks 
are located outside of the leading edges as shown in Figs.(7-ld,2d). On the other 
hand, the entropy production by an O-Type grid is concentrated mainly on the 
small upper region near the tip. Because of this localized high entropy distribution 
the overall flow is not affected so much in comparison with the Fan-Type case except 
for the region of high entropy. And its effect on the shock location seems to be very 
little. 

Again from the comparisons of entropy contours in Figs.(6-2b) and (7-24a,b), 
it can be now seen that the constant entropy lines in Fig.(6-2b) have similar shapes 
to the lower compression parts in Figs.(7-24a,b). This can be expected from the 
fact that the waverider is generated from the elliptic-cone flow. As in Figs.(6-2b), 
in Figs.(7-24a,b) we can observe the vortical layer near the waverider wall where 
the entropy gradient is large. A large entropy gradient can be also seen near the 
sharp tip. The large variation of the entropy there would be a source of error 
and difficulty in the numerical calculation. Above the upper surface we can see 
the entropy change. This is from the shock stand-off at the leading edge. At 
the idealized on-design condition the entropy of the upper region should be the 
freestream value. Recall that any conical stream surface for an elliptic cone in the 
shock layer is composed of stream lines which originate from the same ray on the 
shock surface coming out of an elliptic-cone apex. We know each conical stream 



68 


surface has a constant entropy. Thus, the waverider surface should be a constant 
entropy surface at an idealized on-design condition. But we cannot confirm the 
constant entropy for the waverider surface. The partial reason for this is due to 
the utilization of the finite-volume method (FVM) with its pertinent boundary- 
condition imposition. In the FVM the flow data are not defined at the wall and the 
boundary condition for the pressure is imposed by means of the linear extrapolation 
from the value at the center points of the first grid cells from the wall. A partially 
improved result would be obtained by using a finer grid near the waverider surface. 

7.6 Aerodynamic Forces 

The aerodynamic force coefficients at off-design Mach numbers are shown in 
Figs.(7-25a,b,c) at zero angle of attack for the Type-2? waverider. Both Cl and Cd 
decrease monotonically with increasing The lift-to-drag ratio L/ D decreases 

from 3.52 to 3.27, as M 0 0 increases from 3 to 5. For this Mach number range, 
the variation of aerodynamic forces with respect to freestream Mach number is 
significant. 

Also shown in these figures are the on-design results (M 0 0 = 4) according to 
experiment [6], the full-potential equation [3], and the HSDT approximation for the 
idealized cone-derived waverider (Eq.7-1). At the on-design condition, Cl from the 
Euler equations is about 3.6% lower than the result for the full-potential equations, 
and about 9.6% lower than the result from experiment. The on-design value of 
Cd from the Euler equations is about 13.7% lower than the result for the full 
potential equations, and about 19.5% lower than the result from experiment. The 
experimental result, it should be noted, includes a contribution from viscous effects. 
The L/D ratio according to the Euler calculations are 13.4% greater than for the 
full-potential equations and 7.3% greater than the experimental result. 
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The normal- force coefficient in the y-direction Cm and the axial-force coefficient 
in the i-direction Ca are shown in Figs.(7-25d,e), as a function of a for M = 4. 
The comparison of Cm in Fig.(7-25d) of numerical and experimental data[6] shows 
very good agreement, and Cm varies almost linearly. Figure (7-25e) illustrates Ca 
vs. a. As expected the numerical data are lower than the experimental data, since 
the Ca of viscous flow has an additional contribution due to the skin friction. 

Figures (7-25f,g,h) show aerodynamic quantities Cl, C o, and L/D for the 
elliptic-cone waverider Type-5 as functions of angle of attack a at the design Mach 
number of 4. It is known[22] that the lift/drag ratio of viscous flow for an idealized 
cone waverider, for example, which has infinitesimally thin winglets, approaches the 
value of its corresponding inviscid flow as the cone deflection angle 6 increases. The 
waverider investigated here has 6 = 18.62° for which the difference of lift/drag ratio 
between inviscid and viscous flows is small, that is, the viscous drag is much smaller 
than the wave drag. This indicates that the comparison of the numerical inviscid 
solutions with the experimental viscous results can be justified accordingly. 

Figure (7-25f) shows the lift coefficients Cl as a function of angle of attack at 
the on-design Mach number 4. The agreement of the numerical results obtained 
by solving the complete Euler equations with experimented data is quite good for 
the range of a from —12° and 10°. Approximately, Cl increases linearly with a. 
As a decreases, Cl also decreases and becomes zero at a = —8.5°. For negative 
at, the full-potential result is very close to the experimental data, but the accuracy 
declines as a increases in the positive region. These phenomena can be expected 
from the theoretical restrictions on the potential theory. As a has larger positive 
value, the effective flow deflection angle becomes much larger. This causes the flow 
to be more nonlinear and to increase the rotationality. Accordingly, the homentropic 
approximation which is assumed in the full potential theory becomes worse, as or 
increases. On the other hand, negative angles of attack do not increase the effective 
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flow deflection angle for the range of a used here, considering that at the on-design 
condition the centerline of the waverider (= (5/2) ha$ already positive deflection to 
the freestream direction and the negative angle of attack is compensated by this 
deflection more or less. This is the reason why the agreement is good for negative 
angles of attack. 

The drag coefficients Cd are plotted in Fig.(7-25g) with the same conditions 
as in Fig.(7-25f). The variation resembles a parabolic shape about the a of zero 
lift value. The numerical values are lower than the experimental data, as expected. 
For the positive angles of attack, there is considerable discrepancy between the 
full-potential result and the experimental data, again as for Cl- 

The L/D ratios are depicted in Fig.(7-25h). The maximum value of L/D for 
both the Euler and full-potential calculations occurs at a = —2°, whereas the max- 
imum value for the experimental results occurs at a = 0° (the on-design condition). 
In principle, the L/D for inviscid flow should be greater than that for viscous flow, 
since viscosity increases the drag while affecting the lift only a small amount. The 
Euler results are in accord with this whereas the full-potential results are not, at 
least for a > —2°. Systematic errors in the calculation of Cl and Cd tend to 
compensate when the ratio is taken to obtain L/D. This is especially true for the 
calculation of the reference area since it cancels out entirely for the L / D ratio. 

So far only the Type- 1? waverider has been discussed. To anticipate the L/D 
ratio for another type waverider it would be useful to consider an idealized cone 
waverider which has infinitesimally thin winglets. The L/D for such a waverider 
can be obtained from reference[34] as 

sin <t> 3 4a 2 /(a + 1) 

Hs 1 + ^/na 2 

The Type- .A waverider has smaller 6 and <j> s than Type-B. According to Eq.(7-1) 
these smaller values provide a larger L/D ratio for Type-A. Actually the on-design 
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numerical value of L/D for Type- .A is 6.52, which is much greater than 3.38 for 
Type-B waverider. Those values of L/D due to Eq.(7-1) are slightly less than the 
above numerical results for elliptic-cone waveriders . 



Chapter VIII 

CONCLUDING REMARKS 


A comprehensive study has been carried out for the numerical calculation of 
the inviscid hypersonic flow past a class of elliptic-cone derived waveriders. The 
on-design conditions were Moo = 4 and a = 0°. A variety of off-design condi- 
tions (Moo i=- 4, a ^0°) was studied, all for zero yaw. Thus all the flows were 
symmetric about a central vertical plane. Both the body of the waverider and the 
associated flow fields were conical. An extensive background analysis of the un- 
derlying hypersonic small-disturbance theory was given, and comparative results 
for the on-design flows were presented. The basic flow past the elliptic cone that 
generates the waveriders was also studied as a means for evaluating the flow past 
the waveriders themselves. 

The numerical results for the waverider flow fields for the on-design conditions 
compared well with the HSDT and the numerical elliptic-cone flow fields except 
near the sharp leading edge of the waverider. Here it was found that the bow shock 
stand off from the leading edge rather them being attached as it should according 
to the underlying waverider concept. This was found to be true both for an O-type 
grid and a fan-type grid near the leading edge. The approximate formula use to 
generate the waverider compression surfaces for the A1 and B 1 models was the 
same as was used in previous investigations to study the same waverider shape. 
This approximate formula produces a leading-edge tip angle that is considerably 
larger than it should be to reproduce the proper elliptic-cone flow field. When a 
better approximate formula is used to determine the waverider compression surfaces, 
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giving models A2 and B 2, the correct tip angle is obtained. For model A2, the CFD 
results show that the shock becomes attached at the leading edge, as it should be. 
For model B 2, the shock does not become attached, but the stand-off distance is less. 
For models SI and B 2, the conical shock angle at the tip is about 25°. When the 
angles become this large, the small-angle approximations associated with HSDT 
becomes less accurate. This explains in part some of the remaining differences 
between the HSDT and CFD results, Overall, the agreement of the CFD results 
with what was expected at the on-design conditions was deemed to be good. 

The great importance for the CFD calculations is associated with the off-design 
results. At a = 0° and for M oo > . 4, the bow shock fits more tightly under the 
waverider body. Near the leading-edge of the waverider, a lambda-type shock con- 
figuration appears to develop. This becomes more pronounced as M, oo increases, 
the largest value being M oo = 10 in this study. This is a new effect that was not 
known from previous experimental studies. The lambda-shock configuration occurs 
for other off-design conditions also as part of the adjustment of the main conical 
bow shock to the local leading-edge conditions. 

For Moo < 4 and a = 0°, the bow shock stands off from the body, the more 
so as Moo decreases. There is then a flow in the gap between the shock and the 
leading edge as the flow adjusts to the higher pressure under the body to the lower 
pressure on top. 

When the angle of attack is positive, a > 0°, the shock tends to fit tighter to the 
lower compression surface, and there is an expansion region over the upper surface. 
A very weak shock develops on the upper surface. This is needed to deflect the 
flow that has expanded over the leading edge back parallel to the vertical symmetry 
plane. For negative angles of attack, the shock below the body weakens and a bow 
shock develops over the upper surface of the body. 



74 


The forces on the waveriders could be calculated from the pressure distributions 
on the surfaces. Of course, the viscous contributions could not be calculated from 
the Euler-code results. The results were thus not directly comparable with those of 
experiments, but the trends in variation with angle of attack were very similar. The 
orientation for zero-lift agreed with experiment, and the interpretation in terms of 
the flow-field properties is a unique capability of the CFD analysis. 

For the conical waverider under investigation the Lawrence Code was found to 
be an effective means of producing the numerical results. Care must be taken to 
produce an effective grid structure especially near the sharp leading edges. Unlike 
some other researches for supersonic flows with a sharp leading edge, which ap- 
peared recently, we adopted steady equations instead of unsteady equations for the 
computation. This can save the computer storage and CPU time evidently. Due to 
this saving we could solve waverider flows for a wide range of off-design conditions. 

There is a substantial amount of research that remains for the future. It would 
be useful to calculate the effects of yaw, or sideslip. More important would be 
the calculation of viscous and heat-transfer effects. Attempts associated with the 
present investigation were not successful. In addition, it would also be desirable to 
incorporate real-gas effects. There are other classes of waverider shapes that are 
also of interest. Comprehensive studies relating to their on-design and off-design 
flow fields would be very desirable. 
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Appendix A 

APPROXIMATE ANALYTICAL SOLUTIONS 


The waverider investigated in this study is derived based on the inviscid flow 
field over an elliptic cone. To check the validity of the numerical results of the flow 

over the waverider, it may be necessary to calculate approximate analytical solutions 

/ 

to the inviscid flow around an elliptic cone. The elliptic cone can be considered as a 
perturbed circular cone. Rasmussen et al. [21, 22, 35] obtained approximate solutions 
to the flows around circular and elliptic cones with a small angle of attack by 
means of perturbation theory. They are expressed by linear combination of the 
small perturbed terms of the eccentricity and angle of attack in the framework 
of the Hypersonic Small Disturbance Theory (HSDT). Since the purpose of this 
introduction of the approximate solution is to check the validity of the numerical 
solutions as mentioned before, presented are only the elliptic cone solutions without 
incidence that are corresponding to the flow at the on-design condition of an elliptic 
cone waverider. It will be also useful for understanding the physics of the waverider 
flow to study analytical solutions of an elliptic cone flow. 

The first section of this chapter presents the assumed solutions by outer ex- 
pansion forms. Then unperturbed and perturbed solutions are sought separately 
and they are combined to get the desired analytical approximate solutions. From 
the consideration of streamlines the lower compression part of the waverider con- 
figuration is determined. 
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A.l Perturbation Series Expansions 

For slightly perturbed elliptic cone without incidence, the solutions can be 
approximated by outer expansions as 

u(0, 4> ) = u 0 (9 0 ) + e 2 u 2 (0 o ) cos 2 <f> + 0(e ]), 
v(d , <f>) = v 0 (6 0 ) + £ 2 ^ 2 ( 00 ) cos 2 <f> + 0(e\), 
w(6, <f>) = e 2 w 2 (8 Q ) sin2(f) + 0(e 2 ), 

( A - la, b, c, d, e, /) 

p(9,4>) =p 0 (^o)[l + e 2 P 2 (9o)cos 2<j>] + 0(e\), 
p(9, 4) = po(&o) [1 + e 2 R 2 {0Q) cos 24 } + 0(e \), 
s(6,4) = s 0 + c v e 2 S 2 (0 o )cos2(f> + 0(e 2 ), 

where the subscript 0 denotes an unperturbed basic cone quantity and the subscript 
2 denotes a perturbed quantity. The independent variable 6 0 is defined in terms of 
9 and <f> in a later section. Here P 2 , R 2 , and S 2 axe dimensionless but all the rest 
flow variables axe dimensional. 


A. 2 Unperturbed Basic Circular Cone Solution 

The gas dynamics equation of inviscid axisymmetric cone flow is called Taylor- 
Maccoll equation which can be written as 


( 


2 > 

i-i 

CL 1 . 


cPu 0 duo . fo v o\ n 

_ + cot «_ + 2- ? j Uo = °. 


(A -2) 


The Eq.(A-2) is a highly nonlinear ODE of second order, since both v$ and a 2 are 
functions of uq. The irrotational condition which is obtained from the r- momentum 
equation is given by Vo = duo/d9, and the homenergic condition determines the 
a — uq relation, that is 


a 2 = (7 - 1) 




(A- 3) 


This flow is homentropic for the region between the shock and the cone wall, since 
the entropy jump at the points on the shock is the same. 
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Pottsepp[36] made an approximate relation to Eq.(A-2) which results in a linear 
equation without Uq / a 2 terms. Rasmussen[37j also obtained approximate solutions 
by means of the HSDT. This theory is to assume the limit that 6 — * 0, M 0 0 — ► oo 
, while Ks = Moo sin <5 is held constant. The brief procedure will be introduced in 
the following. First we can see that cos 9 is a solution to the linearized equation. 
This suggests an idea about the Taylor-Maccoll’s solution in the form of 17(0) cos 9 
where U{9 ) is an unknown function of 9. After substitution of the assumed solution 
and neglecting the higher order terms than 0 2 , the desired solutions are obtained 
as 


«»(«) ■ « 2 l-c 2 ^ 

T£- 1 " T ~ 


(A — 4a, b ) 


v a 


where e = p oo/ p(/3 ) is the density ratio across the shock. The gas dynamics equation 
is the second order, but there are three boundary conditions; two at the shock and 
one at the wall. The second order differential equation needs the first two conditions. 
The third one is from the impermeable conditon at the wall. This is necessary for 
the determination of shock location which is another unknown variable. If we apply 
the last restriction to the u-equation, it gives 


where 


P _ 1 

6 y/T=l' 


(.4 - 5) 


(7 - 1 )Kj + 2 

(7 + mi 


( 7 - mi + 2 

(7 + !)*!+ 2’ 


K, = M„I3, 


Ks s M^S. (4 - 6, 6, c) 


After combining the above three relations, the hypersonic similarity shock formula 
can be obtained as 



(4 - 7) 
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Eq.(A-2) was derived by Chemyi[38] and also by Rasmussen by means of the HSDT. 
This represents a very good agreement with the exact formula especially in the 
hypersonic limit. The velocity components uq,vq can be obtained as 


««(»)_, *V 2 ., n A, 

"kT - 1 “ ~2W + fe( F )] ’ 
v~ - 9(1 ¥>■ 


(A — 8a, b ) 


For temperature, we use the homenergic relation like 


T 0 (6) h(8) 


= 1 + 


7-1 


-To o' /loo ' ■ 2 

With the aid of Eq.(A-3) the Eq.(A-9) becomes 

T 0 (e) 


MU 1- 


m 

VI ■ 


(A -9) 


= l + 2 ^K][2 - 6 l + ln(^)}. 


6 2 


(A -10) 


For pressure, homentropic condition can be used to relate the pressure to enthalpy 
which is a linear function of temperature. That is 

'To(d) To(0Y 


Po(fl) = Po(^) + 7 1 


(A -11) 


Poo Poo 7 1 £ L Too Too 

The first term on the right hand side can be obtained through the oblique shock 
relation and the temperature ratio is given by Eq.(A-lO). Thus, we can have 

Po(0) 


Poo 


= 1 + 5*? {i + ;H-£ + '"(£)]}• M-12) 


A. 3 Perturbed Solutions 

The approximate solutions for slightly perturbed cone flows by Rasmussen[22] 
deal with the inclined circular, elliptic, ternary, and quartic cones parametrically. A 
perturbed body can be obtained by identifying a single term out of a Fourier-series 
expansion terms. Here we consider an elliptic cone flows without angle of attack. 
It is well known that the curved shock is a source of rotationality. The elliptic 
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cone has a nonuniform transverse shock curvature and thus the flow is no longer 
irrotational even in the inviscid region unlike the axisymmetric cone flow. The 
analytical solutions about slightly perturbed body can play an important role in 
the studying the physics of the waverider flow where we don’t know exact solutions. 

Ellitic Cone Body and Shock Shape 

A conical body can be expressed by F (9, <f) = 0 as well as its corresponding bow 
shock in spherical coordinate system. For a specific body geometry there must be 
a constraint between 9 and <f> such as 9 = 9 c (<j>). A conical body may be expressed 
by 

F = B - 9 c (<f>). ( A - 13) 

An elliptic cone cross section is given by 

9 c (<f>) = S[l - e 2 cos 2<f> + O(el)}, (A - 14) 

where |e 2 | < 1. In the same line with the body shape, the shock shape may be 
assumed to have the form 9 — 9 S (4>) 

9 3 (<t>) = 8\a - e 2 <72 cos 2 <f> + 0(e 2 )], ( A - 15) 

where g 2 is a shock eccentricity factor which must be determined as a part of 
solutions. 



85 


Boundary Conditions for Body and Shock 

The wall boundary has only tangential velocity with its normal component 
vanishing 

V ■ n = 0, (A - 16) 

where n is an outward unit normal vector. The shock jump conditions can be 
imposed by identifying n at the shock surface. That is, 

n s = eg - e 2 sin (A - 17) 
sinp 

The mass conservation and tangential momentum conservation become the two 
shock boundary conditions. 

V 4 • n 4 = — (Ko ■ n 3 ), (A -18) 

P* 

V s x n 4 = Voo x h a . ( A - 19) 


Stretched Independent Variable 

Considering that 8 is smaller than 6 for some <j) a stretched independent variable 
8 0 is defined by 


_ 8 0 - 6 _ z - 1 _ 8 - 8 c (<j)) 

P-s = 7 ^T = s,w-e c w 

where z = 8 q j8. The three independent variables have the ranges of 


(A - 20) 


8 < 8q < /?, 1 < z < o', 0 < rjo < 1. {A — 21a, b, c) 


On the basis of these independent varibles, the following dependent variables are 
introduced as 

U{8o) = V(z) = U”{t 1 o). (A -22) 


Solution 
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After some algebra and approximation based on the HSDT, the mass conser- 
vation equation can be expressed as 

elu"(e o) + e 0 u'(e Q ) - 4 u{e 0 ) = -2 h 2 (8 0 \ (a - 23 a) 

where 


Ho; 7(7-1)^ 


1 - 


el-s 

p * -S' 2 


(A - 23 b) 


Note that the differential equation is equidimensional and thus the homogeneous 
solution can be obtained. The particular solution can also be obtained by means of 
the method of variation of parameters. The solution in terms of z can be obtained 


as 


CT(*) = Efc(z) + ££(*), 

where its homogeneous and particular solutions are given by 


and 


Up ( z ) J 9i 
W« “ 2a 3 



Z lzl + iflzir 1 _ (ilai)3/2| 
a 2 -l 3 j 2 1 '■<r*-v 1 


z 2 (cos 1 J — 


cos 


-1 1 


i) 


Va 2 — 1 


)' 


where 


J. s 4(6) 


1 + + Ina 2 ) 


4(0) 1 + ^|(2 -1/a 2 )' 


(A - 24) 


(A - 25) 


{A - 26) 


(A - 27) 


The eccentricity factor <72 will be determined at the end of this subsection. Thus 

the perturbed solutions can be obtained as 
u 2 (8q) = U(6 0 ) + uo(#o)02(0o) 5 


v 2 (0q) = U'(8 0 ) + v' 0 (6o)Q 2 (Oo), 

^2(^0) — — — r [2t;o(0o)02(0o) + H 2 (d 0 ) — 2u 2 (6q)\, 
sm uq 

P2(0o) = -7[ u o(0o)tt2(0o) + vo(0o)v 2 (eo)]/al - S 2 (/?)/(7 — 1), 

Ri(0o) = -[u O (0O )u 2 (0 O ) + U 0 (^o) u 2(^o)]/ a 0 - S 2 (P)/{~f - 1), 

(A — 28a, b, c,d, e) 
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where 


ts ta \ — (^0 “ <5)(1 ~ 9l) c c t o\ — /w n* _2 ru 

(J2{('o) = <■>) 02 (p) = 27am 3> a o = 7 — 

a-1 «o(^)/ a 5o a Po 


_ 7(7 - 1)^1 92 n 2 = & 

a l(P)/ a lo <jZ ' 0 7 Po' 

(A — 29a, b, c ) 


From the boundary condition, we have 


V 2 (6) = 0 . 


(A - 30) 


If Eq.(A-30) is combined with ui equation in Eq.(A-28), then g can be obtained as 


1 <7 4 — 1 o- 4 + 1 J f 3 cos 1 i 2 


52 2(J 3 + <r(y + 1) + -6a 3 [ sjo 2 - 1 


(A -31) 


A .4 Streamline and Waverider Configuration 

The compression undersurface of the waverider is obtained from the differential 
equation for streamlines expressed by 


V x dr = 0. 


(A - 32) 


The streamline equation to the first order becomes 

dr rdd 0 rdodcfi 

u o (0o) vq( 6 o) £ 2 ^ 2 ( 00 ) sin 2<£’ 


(A - 33) 


Within the framework of the HSDT with the assumption, 8 0 — ► 8 , the right side 
equation can be integrated to obtain 


#o -6 _ / tan 0 1/e 
/? — <5 tan 0, ’ 


- 34 ) 


where 


ia 2 (<5) *^2(8') _ 2Jg2 2U(8) 

6 = ~ €2 ~8V^' ~8V^ 


(A -35a, 6) 


In the physical variable, Eq.(A-35) become 


(« ) = [ ^ ) + { [?M]_[«M„(W) ,/ e M _ M ) 
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This is another relation describing the compression undersurface of a waverider 
which has a slightly different form from the Eq.(2-13). 

A. 5 Approximate Solution for Shock Laver 

To get the outer solutions of the flow for the compression lower portion of the 
waverider we need to combine Eqs.(A-l), (A-8), (A-10), and (A-12) with Eq.(A- 
28). If we use the constraint on 9 and <f>, Eq.(A-36), the flow variables on the lower 
waverider wall can be obtained. It should be noted that all the approximate analytic 
solutions are not valid across the entire flow field. Nonetheless, the pressure and 
azimuthal velocity component are valid uniformly including the vortical layer which 
exists very near the elliptic-cone body surface. 


Appendix Q 

NUMERICAL ALGORITHM 


The first section introduces a numerical algorithm along with discretization, 
linearization, and operator splitting procedure. The following two sections provide 
the elements of the algorithm; flux Jacobian matrices with the first order expression 
for flux and flux differences with the second order expression for flux. Finally, the 
construction of a block tridiagonal matrix is explained. 

B.l Algorithm 

To numerically integrate the PNS equations we need to express them in dis- 
cretized forms. In this study Lawrence’s algorithm [28] based on Chakravarthy’s 
[19,20] TVD [17,18] scheme is adapted. 

Finite Difference and Finite Volume Methods 

The discretization is accomplished in the frame work of Finite Volume Method 
(FVM) [32]. This utilizes the conservation equations in the integral form and re- 
places the surface integrals with the sum of the flux multiplied by its corresponding 
surface element of a finite volume. There is another discretizing method; Finite 
Difference Method (FDM)[39] which has been widely used. Both of them are in the 
category of local methods. The main differences lie in the treatment of boundary 
conditions and the definition of metrics. The first control points of FDM are lo- 
cated on the boundary, while those of FVM are off the boundary and thus better 
to handle irregular boundary problems than the former. Considering the waverider 
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geometry has sharp comers, the FVM would be more appropriate than the FDM 
for the study of waveriders. The FVM defines metrics at cell interfaces and the 
conservative variables at primary grid points, but the FDM does not separate them 
and use the same points for both metrics and flow variables. 


Discretized Equation in Finite Volume Framework 

Considering that the concept of the physical conservation laws which govern 
the gas dynamics is based on a finite region, let’s express governing equations in an 
integral form for a finite volume V enclosed by a surface 5. 


/ lt dV + / 5 • hdS = J PdV , 


(B- 1 ) 


where U is a conservative variable column matrix and P is the rate of production 
of U. For a steady flow without any source term, it becomes 


j) H • hdS = 0, 


(B- 2) 


where the dyadic H is defined as 


H = Ei 4* Fj 4" Gfc, 


(B- 3) 


and h = outward unit normal vector. Here it is noticeable that for a freestream 
where the fluxes are constant, the integral equation agrees with the geometric con- 
servation, that is, 

j hdS= j dS = 0. (B - 4) 

This means that the integral form of the governing equations satisfies the freestream 
conservation automatically. However, that is not the case for some earlier upwind 
algorithms like that of Buning et al.[40] and Buning[41], and thus various efforts 
have been exerted to remove this drawback. The discretized form of the integral 
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equation for a finite volume which is bounded by six surfaces as in Fig. (5-8) is given 
by 

6 ^ 

^H-AS q = 0, (B- 5) 


where 


0=1 


AS. = AS5*r, AS, S AS£*„ A5 3 =A^j 
AS t = AS^ AS s sA^,, AS t 3AS"+? L , 


p"+ 3 


(B — 6a, b , c,d, e, /) 


and 


— » * A <* 

AS a = / a * + ^oj + n Q k, a = 1,2,..., 6. 


(B-7) 


Note that the subscripts k and l denote primary grids according to Vinokur’s ter- 
minology [32]. Two kinds of grid points such els primary and secondary grids are 
shown in Fig.(5-8). The Eq.(B-5) can be now expressed as 

{£-« _*•} + (jwj -c.-or 5 - °- < b - s > 


with the definition of each element 


4T = H 


A5j, 


f-n+l/2 

fc+1/2,1 


= H • A S 2 , 


f,n+l/2 

^k, 1+1/2 


= H • A5 3 , 



-H ■ A S 4 , 


pn + l/2 
^*-1/2,/ 


-H • A5 5 , 


^n+l/2 
U k, 1-1/2 


= -H • A5 6 , 

(B — 9 a, b, c, d, e, /) 


and 

H • A5 a = l a E + m a F + n a G, a = 1,2, ...,6. (B — 10) 


Eq.(B-8) is a discrete form of Eq.(3-5) for steady flow. At this stage, the above 
expression discretized based on a finite volume frame work is the same as that for 
the finite difference approach. One important difference lies in that l a ,m Q , and n Q 
are calulated at cell interfaces instead of primary grid points. In chapter V it will be 
explained that A S Q are related to themetrics. As in the previous chapter we can 
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obtain a discretized expression for the PNS equations by dropping the streamwise 
viscous derivative terms as 

{^ n+ ' + {^+! + ( 6,+ i - = °- (s - n) 

where E,F, and G are given in Eq.(3-7). 


Streamwise Flow ( Space Marching ) 

We define a streamwise flux Jacobian matrix A * and an operator <5 n+1 as 


A » n = 

A k,l = 


dE*(AS n+1 ,U n ) 


J k,l 


dU 

«■+•() s ()*+■- ()», 

To implement the space-marching we introduce Vigneron’s[42] technique as 


(. B - 12a, b) 


E n = £*(AS n , U n ) + £ p (AS n , t/ n_1 ) 


(B- 13) 


where 


E* = [pU,puU + (jj)up,pvU + (~-)wp,pwU + (y)wp,(F< +p)U] T , 
E” = ( l-^)p[0,(^),(^),(^),01 T , 


with 


u = min 


1 , 




, u.jwv. 


( B — 14 a, 6) 


(J5 — 15a, 6) 


1 + (7-1)A*?j 

where a is a safe factor for nonlinear effect and U is a contravaxiant velocity com- 
ponent in the {-direction. Now the streamwise flux difference can be expressed 
by 


(B - 16) 
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where the marching steps of A* and E p are lagged by one for the conservative vector 
U only. Thus, the operator 6 n in Eq.(B-16) is defined by 

6* K> = , Uh) - A-(ASt,„ Vtj' ), 

(B — 17a, 6) 

= ^( A ^+ 1 , OtJ) - E”(ASl„ VZJ'). 


From now on, the E* replaces the role of E and the E p is neglected or treated like 
a source term in the space marching scheme. 


Crosswise Flows ( Linearization ) 

To avoid confusion and make the linearization procedure of cross flows clear, 
some notations as in Table 3 are introduced. 


Table 3 

K 

m 

H 

Hi 

Hv 

V 

k r 

F 

Fi 

F v 

C 

k\\ 

G 

Gi 

Gv 


where H = H i + H v and the superscript * is for an inert index by an opreator, which 
means that the index with * does not vary. The H denotes a difference expression of 

A . 

an inviscid and/or a viscous flux. But H is strictly reserved for the flux itself such 
as either F or G. The subscripts „ stand for the inviscid and viscous respectively. 
The general flux difference is defined by a difference operator 6 K as 

HIT 1 = H”-» /2 - h; + J 1/2 , (B - 18) 

where we raised the marching level (n + 1/2) to (n + 1) for the purpose of both 
simplicity and stability. In the above equation the subscript m denotes the cell 
center point and the subscripts m ± 1/2 denote the cell interface points. Before 
starting linearization it is to be noted that 

A — 

1/2 = fn(Um,U m ±i, AS m ±i/2)' 


(B- 19) 
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This means that the flow quantities are defined at the cell centers, while the metrics 
are defined at the cell interfaces where the fluxes are calculated also. Considering 
the elements of in Eq.(B-19) we can linearize it as 


3H n 3H n 

Hit 1 , 3 H! ±1/2 + + .n ±l ' 2 6 n+1 u m±u 

OU m OU m ± i 


L m±l/2 


(B - 20) 


where 


S n+1 U = U n+1 - U n . 


(B- 21 ) 


Depending on the above linearization the general flux difference at the marching 
level (n + 1) can be now obtained in terms of n marching level values 


<„H: +1 = <5„H! + 


* 


dU \ 


(6 K H m ) + 6 k [ d - 


6 n+1 Um, 


where 


c tt n Tjn Tjn 

= n m+ 1/2 “ n m-l/2i 


5H" 






(B - 22) 


(. B - 23) 


d - , - 1/2 _ 

(fl - 24) 

OUm dUm UUm 


a S^Um = Qr -? +1/ ^ n+1 ^ m+1 - ”- 1/8 ^ 1 & m - 1 . (5 - 25) 


a£T m r ^ m+1 at/ m _ a 

It is to be noted that the S K in the Eq.(B-25) operates on the whole portion which 
follows the operator. 


First Form of Algorithm 

To have a desired discretized algorithm all the terms containing 6 n+1 U are 
moved to the left hand side and the other terms to the right hand side of the 
equation after substituting the appropriate expressions for (5 K HJ^ +1 into Eq.(B-ll) 
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with the help of the Table 3. Then the following algorithm in the so called (5-form 
can be obtained, 


'A-+jL {8n r+ k G’ )+ 6 n %+k% 


6 n+ 'Uk.i 


-I *,1 


(B - 26) 


= - {((5 s i*)i/ + 6 h EP + SjjF 11 + 6 c G 7/ }” t , 


where []£ j and {}£ f are for a Jacobian matrix and a flux respectively and 


F 7 = F 7 + F v , G 7 = G 7 + G v , 

F 77 = p.i + F v , G 77 = G 77 + G v . 


(B - 27a, b) 


The superscripts I and II denote the first and the second order expressions for the 
inviscid fluxes respectively. More generally, we introduce the notations , H 7 and 


H 77 , like 


H 7 = either F 7 or G 7 , 


(B — 28a, b) 

H 77 = either F 77 or G 77 . 

It must be noted that these notations are not the general fluxes, but they are the 
expressions for the first and the second order difference representation of them. They 
will be specified in the following sections. If the first order accurate discretization 
is used on the LHS, then block tridiagonal matrices are obtained. If we choose 
the second order expression, then it will become pentadiagonal matrices which will 
require much more efforts to inverse them than the tridiagonal matrices. In view of 
the above fact it would be adequate to use the first order differencing the LHS. On 
the other hand, the RHS must be handled with much care and a lot of efforts are 
to be exerted to imporve the accuracy. For the RHS the second order differencings 
are usually used. 
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Final Form of Algorithm (Approximate Operator Factorization) 

The Eq.(B-26) cannot be used directly. To numerically integrate this it is 
rewritten by approximate operator factorization [43] as 




(B - 29) 


where 

Ws 
\ 0 } = 


A* + i^(Ff + F v ) + F[ + F v ) 


k,t 


A* + + G„) + + G v ) 


(B - 30a) 
(B - 306) 


J k,i 


{RHS} n = — 


( B - 30 c) 


>(A5”r , U? tl ) - A*(AS£ l , ^J 1 )] 

- {E p (AStf\ UZt) - E p (ASl h U?J')} 

At this stage any reasonable finite difference expression can be chosen for H 7 or 
H 77 . The special form used for this study will be presented later. The solving 
procedure of the Eq.(B-29) can be taken as the following four steps. 

1. [a}X = {RHS} n , 


2. Y = 


A*n 

A k,l 


X , 


3. [0}8 n+1 U k ,, = Y, 

4 . =Ki+ 6n+ 'v>.j, 


( B — 31 a, 6, c, d) 


Both of [ a ] and [0] is a block tridiagonal matrces which will be explained in the last 
section. 
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B.2 Flux Jacobian Matrices 

The LHS of the previous algorithm Eq.(B-29) contains some Jacobian matrices 
in [a] and [0]. Their calculating procedure is described in the following subsections. 



where R~ l and R are the left and right eigen vector matrices of the Jacobian 
matrix C respectively. To evaluate R~ l , JZ, and A Roe-averaged values are used. 
The notation “ " ” is introduced to denote the Roe average values. The first order 
expression of the general inviscid flux can be now expressed in terms of sgnC by 

(HfWi/2 = 5 {(ffiWl + (&)«} - \(sgnC) m+U2 AH h (B - 34) 

where 


A() = () m+1 -() m , \C\ = R(A + -A~)R-\ sgnCAHi = \C\AE*. 

(B — 35 a, b , c) 

Based on the Eq.(B-34), the inviscid Jacobian matrices are calculated as 


a( aoLT. 1/a = (^'•CW I/2 ]B j (A5 m+1/2 ,if m+I ), 

d ad*' n = 


(B — 36a, b) 
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with the definition of B< as 


^i(A5 m +i/2,[/ m +i) = 


_ a ( H ,-) m+1/2 


au . 


m+1 


Bi(AS m +i/2, U m ) = 


_ 0(Hj) m+1/2 


dU„ 


( B — 37 a, b ) 


For the case of (Hf) m _j/ 2 m in Eqs(B-36,37) is to be changed into m — 1. 


Viscous Flux Jacobian 

For viscous flux Jacobian matrices they are simply defined as 

B v (As m+1/2 ,u m+1 ) =- ( f;- )m+ -^, 

at/ rn+l 

n / a ? fr \ — d(Hv)m+l/2 

i>m+l/2,Um) = -qq • 


(. B — 38 a, 6) 


Inviscid and Viscous Flux Jacobian 
Considering the definition of H 


H 7 = Hf + H v . 


(. B - 39) 


Jacobian matrices for both the inviscid and viscous terms axe obtained from the 
previous two subsections as 

A / 

agT' 78 = 5 [I± (^"<?) m±1 / 2 ]B i (AS m±1/2 ,£? ra ) T B v (AS m ± 1 /2,U m ), 

(B — 40a, b ) 

where the matrices C, I?,-, and B v axe given in reference[28]. The above two equa- 
tions can determine all the terms on the LHS of Eq.(B-29) except for A* which is 
given already in Eq.(B-12). 
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B.3 Flux Differences 

The RHS of the final algorithm of Eq.(B-29) contains the flux differences of the 
second order in F and G as can be seen in Eq.(B-30c). The general flux difference 
is divided into the inviscid and viscous portions as 

6 K (k n )m m S.( Hl') m + «,(£„)„• (B - 41) 

The inviscid flux difference is defined as 

*«( Hf')„ = (H?') m +i/2 - (H.Vi/1- (B - 42) 

The second order expression for Hi [28]is obtained by adding higher order correction 
terms to the first order expression in Eq.(B-34) like 

(H")m±l/ 2 - AF3) m±1/2 

.. (B — 43) 

+(^){AF+ - AF 2 } m±1/2 . 

Depending on the value of the accuracy parameter 4> various schemes can be ob- 
tained. For example, if <f> = 1, it becomes an upwind scheme. The notations ~ and * 
denote the limited flux differences through the MINMOD limiter which is defined 
by 

MINMOD(x,y ) = sgn(x) • max{0,min(\x\, fisgn(x)y)}, ( B — 44) 

where /3 is a compression parameter. The MINMOD picks up a minimum absolute 
value (modulus), if the two arguments are of the same sign. On the other hand, 
it assigns zero, if they are of opposite sign. Therefore, it plays a role of limiting 
the fluxes, whenever it meets very steep varations. But for mild variation it returns 
un-limited fluxes, since is large compared with 1 in most case. At maximum, min- 
imum, or inflection points where the sign of the neighboring fluxes are changed, the 
limiter returns zero. As a result, if we use MINMOD based on the Total Variation 
Diminishing (TVD) condition, it will not introduce new maxima or minima and 
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thus it removes glitches and expansion shocks. This is somewhat parallel with the 
Piecewise Parabolic Method (PPM) by Collela and Woodward[44]. As /5 becomes 
larger, we can get a sharper distribution near discontinuity. However, to satisfy 
the TVD condition, its value is restricted. For instance, in one-dimensional scalar 
nonlinear equation or systems of linear equations, 

l</?<l=i. (B — 45) 

Even though the equations used here are multi-dimensional non-linear systems of 
equations, it is known that this can provide quite good guideline for them also. 
The viscous flux difference which can be simply determined by 

<5 < c(B r t))m = (H v )m+ 1/2 — {B — 46) 

B.4 Block Tridiagonal Matrix 
Matrix Construction 

In order to integrate the PNS equations by the algorithm given by Eq.(B-29) 
the Eqs.(B-40) and (B-38) must be inserted in Eq.(B-30). In doing so the first step 
of Eq.(B-31) is important, since it contains a block tridiagonal matrix[45] on the 
left hand side. The third step is the same as the first in the point of numerical 
calculation. The rest two steps are not difficult. Thus, the solving procudure of 
block tridiagonal matrix will be presented in the following. The Eq.(B-31a) can be 
rewritten as 




(B - 47) 
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where 




d(F J ) t -i/ 2 ,j 

dUk- 1,1 ’ 


[ck2]jfc,/ = 


_ d(r) k+1/2tl 5(F0 t -i /2l , 


dU k ,i 


+ (B-48a,b,c) 


d(F f )t +1 / 2 ,I 

= -sscr; ' 


Here the [ai]t,i, [02]*, h erid [ 03 ]*,! are 5x5 square matrices. {XJr-ij, { A" } t, 
{*}*+■,!, and {RHS}k,i are 1x5 column matrices. Thus the coefficient matrices, 
[a*], compose a block tridiagonal matrix. Eq.(B-47) is solved along the rj -direction 
( / = const ), while Eq.(B-31c) is solved along the (^-direction ( k = const ). Thus the 


calculating procedure changes the directions. For this reason it is called Alternating 
Directional Implicit (ADI) method. 


Matrix Solver 

A system of equations composed of a scalar tridiagonal matrix can be solved by 
means of Thomas algorithm. For a block tridiagonal matrix, we can also apply the 
same algorithm combined with the inverse matrix solver. The solution procedure of 
a system of equation with tridiagonal matrix is described in brief. First, we eliminate 
the sub-diagonal matrices below the diagonal matrices by using the upper row 
matrices and then the present row matrices are multiplied by the inversed diagonal 
matrix element. Thus we can construct upper triangular matrix whose elements 
are 5x5 block matrices. To obtain the inverse matrix we use Crout’s algorithm[46] 
which decompose any matrix into lower and upper triangular matrices. The final 
step is to take the backward sweep to get the desired block column matrices which 
contain the dependent (conservative) variables. 



Appendix C 

GEOMETRY OF LEADING EDGE 


It is desired to determine the angle A between the upper planar surface and 
the lower compression surface. The outward pointing unit normal of compression 
surface is determined by n = — VF/|VF|, where F(0, <j>) = 6 — 6(<j>), where 6(4>) is 
determined by Eq.(2-6). For small polar angles, we have 


n = 


z B _ 
e 6 ~ g ] 


1 + fliiV 

1 + ^0 d<t> 


) 


1/2 ' 


(C- 1 ) 


Since the outward normal to the freestream planar surface is e^, we can determine 
the leading edge tip angle A (where 6 = 0) by cos(7r — A) = n • e^. It can then be 
determined that 


tan A = P (-jz 


£2^2(0) sin 2(f> s 


d6 J 


(C-2) 


VooMi - 7*) 

The azimuthal velocity at the shock can be determined from the shock boundary 
conditions [21,22] : 

M/?) _ 2 52 (c _ 


SV„ a* 

Thus the leading-edge tip angle is determined by 

2e2<72 sin 2 cf> s 


tan A = 


(C-4) 


a(o 2 — 1) 

The tip angle according to the approximate formula (2-9) is determined by a differ- 
ent result. It is 


(tan A) 


eo 


approx 2 (cr — 1 ) 


sin2<^> 3 . 


(C — 5) 
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The tip angles for models A1 and B 1 are determined by the approximate 
formula (C-5), whereas the tip angles for models for A2 and B 2 are determined by 
the correct formula (C-4). The numerical values are shown in the following table: 


Table (C-l) 


Model 

A 

Al 

10.45° 

A2 

1.44° 

B 1 

9.45° 

B2 

4.12° 


The models A1 and B 1 with the compression surface described by the approxi- 
mate formula (2-9) have considerably thicker angles than the corresponding correct 
waveriders would have. 
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Fig.(l-l) Model of Aero-space Plane 
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Fig.(2-1) Spherical Polar Coordinates 
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Fig.(2-2) Construction of Elliptic-Cone Waverider 
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y/ztand 

Fig. (2-4) Waverider Type— A 
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(x, y, z) = Cartesian 
(r, 0, <j)) = Spherical 
{%, T|, Q = Body-fitted 


Fig.(3-1) Coordinate Systems 
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Fig.(3-3) Boundary Conditions 
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Fig. (3-4) Reflection Condition 







114 



C 


max 


( Physical Domain ) 

c 



( Computational Domain ) 


Fig.(5-1) Elliptic Grid Generation 



Fig. (5-2) O-Tvpe Grid 
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Fig. (5-3) Magnified 0-Type Grid 



Fig. (5-4) Fan-Type Grid 
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Fig. (5-5) Magnified Fan-Type Grid 




Fig. (5-6) Adaptive Grid 
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Fig. (5-7) Magnified Adaptive Grid 
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• = Primary Grid Points 
X = Secondary Grid Points 


Fig.(5-8) Finite-Volume Element 
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Fig. (6-la) Grid for Elliptic-Cone Flow ( for WR Type-R ) 
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Fig. (6- If) Wall Pressure Comparison 
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y/xtand 


Fig. (6-lg) Comparison of Pressure near Lower Symmetry Plane 
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Fig. (6-lh) Pressure Contours 
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Fig. (6-li) Mach Number Contours 





Fig. (6-2a) Entropy Contours ( <5 = 12° Elliptic Cone ) 






Fig. (6-2b) Entropy Contours ( <5 = 18.62° Elliptic Cone ) 
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Fig. (6-3a) Shock Comparison ( 8 = 12° Elliptic Cone ) 
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Fig. (7-la) Grid for Type- .410 Waverider 


( On-Design, WR Type-AlO ) 
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Fig. (7-lc) Magnified Cross-Plane Velocity near Tip 


( On-Design, WR Type-AlO ) 
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Fig. (7-ld) Shock Location Comparison 
( On-Design, WR Type-AlO ) 
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Fig. (7-lg) Pressure near Lower Symmetry Plane 


( On-Design, WR Type-AlO ) 


Fig. (7-lh) Pressure Contours 
( On-Design, WR Type-AlO ) 



Fig. (7-li) Mach Number Contours 


( On-Design, WR Type-AlO ) 
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Fig. (7-2a) Grid for Type-^42i ? Waverider 



Fig. (7-2a') Magnified Fan-grid for Tip Region 
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Fig. (7-2c) Magnified Cross-Plane Velocity near Tip 
( On-Design, WR Type-A2F ) 
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Fig. (7-2d) Shock Location Comparison 


( On-Design, WR Type-A2F ) 
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Fig. (7-2e) Azimuthal Velocity Component 
( On-Design, WR Type-A2F ) 
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Fig. (7-2f) Wall Pressure Distribution 


( On-Design, WR Type-A2F ) 







FlS ’ (7-3a) Grid for Type-BlO Waverider 
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Fig. (7-3c) Magnified Cross-Plane Velocity near Tip 
( On-Design, WR Type-BlO ) 
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Fig. (7-3d) Shock Location Comparison 
( On-Design, WR Type-BlO ) 
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Fig. (7-3e) Azimuthal Velocity Component 
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Fig. (7-3g) Pressure near Lower Symmetry Plane 


( On-Design, WR Type-BlO ) 



Fig. (7-3h) Pressure Contours 
( On-Design, WR Type-BlO ) 




Fig. (7-3i) Mach Number Contours 
( On-Design, WR Type-BlO ) 
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Fig. (7-3k) Mesh Size Effect ( Shock ) 
( On-Design, WR Type-BlO ) 
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z/xtand 


Fig. (7-31) Mesh Size Effect ( Leading Edge ) 
( On-Design, WR Type-BlO ) 




Fig. (7-4a) Grid for Type-B2F Waverider 
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Fig. (7-4b) Cross-Plane Velocity Distribution 
( On-Design, WR Type-B2F ) 
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Fig. (7-4c) Magnified Cross-Plane Velocity near Tip 
( On-Design, WR Type-B2F ) 
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Fig. (7-4d) Shock Location Comparison 
( On-Design, WR Type-B2F ) 
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Fig. (7-4e) Azimuthal Velocity Component 
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Waverider Type-AiO 





Fig. (7-5a) u : Radial Velocity Component 
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Waverider Type-AiO 
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Fig. (7-5b) v : Polar Velocity Component 
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Fig. (7-5c) w : Azimuthal Velocity Component 
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Waverider Type-AlO 



Fig. (7-5d) p : Pressure Distribution 
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Waverider Type-AiO 



Fig. (7-5e) p : Density Distribution 
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Waverider Type-AlO 
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Fig. (7-5f) T : Temperature Distribution 




oo 


177 


Waverider Type-BiO 
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Fig. (7-6a) u : Radial Velocity Component 
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Fig. (7-6b) v : Polar Velocity Component 
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Waverider Type-BlO 
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Fig. (7-6c) w : Azimuthal Velocity Component 
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Fig. (7-6d) p : Pressure Distribution 





Fig. (7-6e) p : Density Distribution 
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Waverider Type-BiO 
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Fig. (7-6f) T : Temperature Distribution 
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Fig. (7-7b) Wall Pressure Distribution 


(Moo = 3, a = 0°) 
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Fig. (7-7e) Pressure Contours 
(M 0 o =3, a = 0°) 
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Fig. (7-7f) Mach Number Contours 
(Moo = 3, a = 0°) 
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Fig. (7-8b) Magnified Cross-Plane Velocity near Tip 
(Moo = 4.5, a = 0°) 
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Fig. (7-8c) Wall Pressure Distribution 


(Moo = 4.5, a = 0°) 
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Fig. (7-9b) Magnified Cross-Plane Velocity near Tip 
(Moo = 5, * = 0°) 



Fig. (7-9c) Total Pressure Contours 
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Fig. (7-9d) Wall Pressure Distribution 
(Mqo = 5, a = 0°) 
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Fig. (7-9g) Pressure Contours 
(Moo = 5, a = 0°) 





Fig. (7-9h) Mach Number Contours 
(M* = 5, a = 0°) ■ 
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Fig. (7-10b) Total Pressure Contours 
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Fig. (7-10c) Schematic Diagram of Lambda Shock 
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Fig. (7-10d) Wall Pressure Distribution 
(Moo = 10, or = 0°) 
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Fig. (7-1 Of) Azimuthal Velocity Component 

(Moo = 10, a = 0°) 
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Fig. (7-1 lb) Wall Pressure 
£ 4, a = 0°) 
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Fig. (7-12a) Cross-Plane Velocity Distribution 
(Moo =4, a = +3°) 
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Fig. (7-12b) Magnified Cross-Plane Velocity near Tip 

(M co = 4, <* = +3°) 
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Fig. (7-12c) Wall Pressure Distribution 
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Fig. (7- 13a) Cross-Plane Velocity Distribution 
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Fig. (7-13b) Magnified Cross-Plane Velocity near Tip 

a = + 10 °) 
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Fig. (7-13c) Wall Pressure Distribution 
(Moo = 4, a = +10°) 




219 



.56 



yjx tan <5 


Fig. (7-13e) Pressure near Upper Symmetry Plane 
{Moo = 4, a = +10°) 
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Fig. (7-14b) Magnified Cross-Plane Velocity near Tip 
(Af«, = 4, a = -2°) 
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Fig. (7- 14c) Wall Pressure Distribution 
(M=c =4, a = -2°) 
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(Afoo = 4, a = -2°) 
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Fig. (7-15b) Magnified Cross-Plane Velocity near Tip 
(Af«, = 4, a = -4°) 
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Fig. (7-15c) Wall Pressure Distribution 


(Moo = 4, a = -4°) 
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Fig. (7-15e) Pressure near Upper Symmetry Plane 
(Moo = 4, 


a = -4°) 
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Fig. (7-16b) Magnified. Cross-Plane Velocity near Tip 
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Fig. (7-16c) Total Pressure Contours 
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Fig. (7-16d) Wall Pressure Distribution 
(Moo =4, a = -8°) 


i 

1.50 











■ : .'////// ////yy >“ \ 

• : / • /: ^^//////.^ 

: • :' : - 




1 

!S®Mj i it f. "‘•iiniiiii 


ill 1 ‘ \ \ \ \ * ■ 


i 1 , * i i * , / 

i 1 , 1 i i y 

* i 1 i i / 


. 1 , 1 . 

* 1 , 1 < 

■ i 


• , * « * * * * , i • 

i , 1 * i i • ' • * , i 


‘ , 1 1 I 1 ' 1 I 1 I 1 1 * , v 

‘ , ‘ , 1 . l . V 

1 l I , 1 t I I l I I 1 / 

• ' ‘ • , 1 - ‘ / 


1 i i i i y 

• i i i \ / 

1 i t j/ 

1 1 1 i/ 


Fig. (7-17a) Cross-Plane Velocity Distribution 
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Fig. (7-17b) Magnified Cross-Plane Velocity near Tip 


(M a o = 4, 


a = -12°) 
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Fig. (7-17c) Wall Pressure Distribution 
(Moo = 4, oc = —12°) 
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Fig. (7-17e) Pressure near Upper Symmetry Plane 
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Fig. (7- 18a) Pressure near Lower Symmetry Plane 
(Afoo = 4, a jk 0°) 
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Fig. (7-18b) Wall Pressure 
(Moo = 4, a ^ 0°) 




Fig. (7- 19a) Cross-Plane Velocity Distribution 
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Fig. (7- 19c) Wall Pressure Distribution 
(Moo = 5, a = +4°) 
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Fig. (7-19e) Pressure near Upper Symmetry Plane 
(Moo = 5, a = +4°) 
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Fig. (7-20b) Magnified Cross-Plane Velocity near Tip 
(Moo = 5, a = -4°) 
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Fig. (7-20e) Pressure near Upper Symmetry Plane 


(Mcc = 5, 


a = -4°) 
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Fig. (7-2 lb) Magnified Cross-Plane Velocity near Tip 
(Moo =3, a = +4°) 
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Fig. (7-21e) Pressure near Upper Symmetry Plane 


(Moo — 3, oc = +4°) 
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Fig. (7-22b) Magnified Cross-Plane Velocity near Tip 


(Moo = 3, 


a = -4°) 
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(Moo = 3, a = -4°) 
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Fig. (7-22e) Pressure near Upper Symmetry Plane 


(Moo = 3, 


a = -4°) 







273 



y/xtand 

Fig. (7-23a) Pressure near Lower Symmetry Plane 
(Moo £ 4, c* # 0°) 


Lower Surface 



“ Mach=3.0, Alpha = -4 
"Mach=3.0, Alpha= + 4 
"“Jtfach=4.0, Alpha=0 
~Mach=5.0, Alpha=-4 
~ Mach=5.0, Alpha^Hf 





Fig. (7-24a) Entropy Contours for Waverider Type-BIF 1 


ENTOOPV 

Uauarldar t Type-BIF 





Fig. (7-24b) Entropy Contours for Waverider Type-510 





2 


Waverider Type-BlO 


0.90 f 



Fig. (7-25a) Lift Coefficient vs. M 






279 


Waverider Type-BlO 



\ 

i 

Euler 



Fig. (7-25c) Lift/Drag Ratio vs. Moo 



Waverider Type-BlO 



Fig. (7-25d) Normal-force Coefficient vs. a 
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Fig. (7-25e) Axial-force Coefficient vs. a 
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Fig. (7-25h) Lift/Drag Ratio vs. a 







